Calculation of binding energies and masses of quarkonia in analytic QCD models 
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We extract quark masses m q (q = b, c) from the evaluation of the masses of quarkonia T(1S) and 
J/if)(lS), performed in two analytic QCD models, and in perturbative QCD in two renormalization 
schemes. In analytic QCD the running coupling has no unphysical singularities in the low-momentum 
regime. We apply the Minimal Analytic (MA) model of Shirkov, Solovtsov, et al.; and the 2-delta 
analytic model (2danQCD). The latter, in contrast to MA, at higher energies basically coincides with 
the perturbative QCD (in the same scheme). We use the renormalon-free mass rn q as input. The 
separation of the soft and ultrasoft parts of the binding energy E qq is performed by the requirement 
of the cancellation of the leading infrared renormalon. The analysis in 2danQCD indicates that the 
low-momentum ultrasoft regime is important for the extraction of the masses m q , especially nib. 
2danQCD gives us clues on how to estimate the influence of the ultrasoft sector on m q in general. 
These effects lead to relatively large values nib ~ 4.35 GeV in 2danQCD. In perturbative QCD in 
MS scheme these effects are even stronger and give larger uncertainties in nib. The MA model gives 
small ultrasoft effects and the extracted values of mj, agree with those in most of the literature 
{nib ~ 4.2 GeV). The extracted values of fn c in all four models are about 1.26-1.27 GeV and agree 
well with those in the literature. 

PACS numbers: 12.38. Cy, 12.38.Aw,12.40.Vv 



Most of the calculations of the masses of heavy quarkonia are based on perturbative expansions. These expansions 
come from the knowledge of the static quark-antiquark potential V(r) which has been calculated up to N 2 LO (~ ce^/r) 
in Refs. [H[2]. At N 3 LO (~ af/r) level, ultrasoft gluons contribute 0H] to V(r) and the calculation of these terms has 
been completed in Refs. 016]. The expansion coefficients for the mass of the (heavy) qq quarkonium vector (S = 1) 
or scalar (S = 0) ground state (i.e., with: n = 1 and I = 0) are given, for example, in Ref. [7J, for the terms including 
~ m q a^. The latter term is completely known because the parameter 0,3 in the corresponding aforementioned ~ 
contribution to the static potential is now known, Refs. For a review of the topic, see, e.g. Ref. [5]. 

Most of the radiative contributions to the quarkonium mass expansion are from the so called soft sector of gluon 
momenta, Q ~ m q a s } In the case of bb and cc, this is around 2 GeV and 1 GeV, respectively. At Q « 2 GeV 
scales, the perturbative QCD (pQCD) coupling a s (Q 2 ) is marginally reliable; at Q« 1 GeV it is unreliable. This is 
so because the pQCD coupling a s (Q 2 ) suffers from unphysical (Landau) singularities at low spacelike q 2 (= — Q 2 ), 
i.e., at < Q 2 < A 2 where A 2 ~ 10 _1 GeV 2 ; for Q 2 ~ 1 GeV 2 it is dangerously close to these singularities, thus 
unreliable. 

The aforementioned Landau singularities are not physical because they do not possess the analytic properties of the 
spacelike observables T>(Q 2 ) (such as Adler function), the latter properties being dictated by the general properties 
of quantum field theories (9] [10] including causality. Namely, T>(Q 2 ) must be an analytic function in the entire 
Q 2 -complex plane, with the exception of the timelike semiaxis Q 2 < —M 2 hr , where M t h r ~ 10 -1 GeV is a particle 
production threshold. Specifically, the evaluation of spacelike quantities T>(Q 2 ) in pQCD, as a (truncated) power 
series of the perturbative running coupling (couplant) a v t(nQ 2 ) = a s (KQ 2 ) fit (with k ~ 1), does not respect these 
important analytic properties of T>(Q 2 ). These problems, within the context of QCD, were first addressed in the 
seminal works of Shirkov, Solovtsov et al. [TlT - TH] . There, the perturbative running coupling a p t(Q 2 ) = a s (Q 2 )/ir was 

made an analytic (in the aforementioned sense) function of Q 2 , a p t(Q 2 ) i-> A^ A \Q 2 ), by keeping in the dispersion 

relation the discontinuity function Pi\cr) = Im a pt (Q 2 = —a — ie) unchanged on the entire negative axis in the 
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1 The momentum transfer q of the gluon is spacelikc, i.e., q 1 = —Q 2 < 0, because the scattering of the quark and antiquark in the 
quarkonium is of the t-channel type. 
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complex Q 2 -plane (i.e., for a > 0), and was set equal to zero on the nonphysical cut < Q 2 < A 2 . For this reason, 
this analytic QCD (anQCD) model will be referred to here as Minimal Analytic (MA). 2 

Various other analytic QCD (anQCD) models for Ai(Q 2 ) can be constructed, and have been presented in the 
literature, among them Refs. [T51423J . These models fulfill certain additional constraints at low and/or at high Q 2 . For 
further literature on various analytic QCD models, we refer to review articles Refs. [24-27 . Some newer constructions 
of anQCD models of A\{Q 2 ) include those based on specific classes of f3 functions with nonperturbative contributions 

[28] or without such contributions [29l [30] , and those based on modifications of the discontinuity function pf 3 i-> p\ 
[ pi(u) = Im A\(Q 2 = —a — ie) ] at low (positive) a where p\ is parametrized in a specific manner (with delta 
functions), cf. Refs. [SI [32J. 

The model of Ref. [32], an extension of the model of Ref. [31], will be called 2-delta analytic QCD (2danQCD), 
because its discontinuity function px(cr) is parametrized in the unknown low-cr regime by two delta functions. It has a 
specific attractive feature that in the high-scale regime (\Q 2 \ > A 2 ) it practically coincides with the corresponding (in 
the scheme) pQCD coupling: A\{Q 2 ) — a p t{Q 2 ) ~ (A 2 /Q 2 ) 5 . This makes it possible for the model to be applied with 
Operator Product Expansion (OPE) [33] including terms of dimension D < 8 without the need to modify the ITEP 
School interpretation of the OPE [331 133 ■ The latter basically states that the terms in OPE of higher dimension 
(D > 0) originate from the infrared regime only. 

Once having an anQCD model, i.e., a model for the analytic analog A\(Q 2 ) of the pQCD coupling a p t(Q 2 ) (or, 
equivalently, for the discontinuity function pi(cr)), the analytic analogs A n (Q 2 ) of the higher powers a pt (Q 2 ) n have to 

be constructed from it. In the case of MA/APT (4 MA) (<2 2 ))> the couplings A^iQ 2 ) were constructed in Ref. 136 
for integer n and in Refs. [37-39 for noninteger n (cf. also Ref. [27]). In the general case of analytic A\{Q 2 ), the 
quantities A n (Q 2 ) were constructed in Refs. [20l [21] for integer n, as linear combinations of logarithmic derivatives 
Ak{Q 2 ) oc d^ 1 Ai(Q 2 ) /d(lnQ 2 ) fc ~ 1 (k > n), 3 and were extended to noninteger n in Ref. [41]. 

In this work, we apply two anQCD models, MA/APT of Refs. [TTHTi] . and 2danQCD of Refs. [31 123], to evaluations 
of the perturbation series of the binding energy E q q of heavy quarkonia (T(1S) and J/ip(lS)) and of the quark pole 
mass ffiq , and thus to evaluations of the masses of these quarkonia M qq = 2m q + E qq . As input parameter we use the 
renormalon-free quark mass fn q (= m q (p 2 — fn 2 )) of the corresponding quark q — b,c (also called MS quark mass), 
and the anQCD coupling of the model. Since the quarkonia masses are well measured, we can extract the values of 
fn q . We also perform the same analysis in pQCD in the corresponding renormalization schemes. 

In Sec.[n]we briefly describe the MA/APT model (Subsec. [!lA| and the 2danQCD model (Subsec. pB] ). In Seep! 



we present the procedures of evaluation of the binding energy E qq and of the quark pole mass m q , in terms of the 
mass ffiq and of the couplings. Furthermore, we explain how the cancellation of the leading infrared renormalon in 
the sum 2m q + E q q allows us to separate the ultrasoft from the soft part of the binding energy. The numerical results 



and the extractions of the masses mb and m c are presented in Sec. IV (Subsections IV A and IV B respectively) 



The evaluations are performed in the two aforementioned analytic models MA/APT and 2danQCD, and in pQCD 
in the corresponding two renormalization schemes (MS, and in the scheme of 2danQCD called here the Lambert 
scheme). In Sec. [V] we summarize our results and draw certain conclusions. Appendix [A] contains the expressions of 
the coefficients of the perturbation expansion of E qq , available from the literature. Appcndix[B]contains a renormalon- 
based estimation of the coefficient at a^ t term in the perturbation expansion of m q /ffi q . 



II. TWO ANALYTIC QCD MODELS 

A. Minimal Analytic (MA/APT) model 

As already mentioned in the Introduction, the coupling a pt (Q 2 ) = a s {Q 2 )/iT in pQCD in the usual renormalization 
schemes (such as MS, or the = /3a = • • • = scheme), possesses unphysical (Landau) singularities inside the required 
analyticity regime Q 2 € C\(— oo, 0]. Specifically, it has a cut on the semiaxis (— oo, A 2 ) (where A 2 ~ 10 _1 GeV 2 is the 
"Landau" branching point), thus offending the analyticity requirement on the cut sector Q 2 £ (0, A 2 ). Application of 



2 In the literature, it is also referred to as Analytic Perturbation Theory (APT). We refer to it as MA, or MA/APT. 

3 The relations between such functions A-k's and ^4 n 's, allowing a recurrent construction of An, for integer n within the context of the 
MA/APT model, were given also in Refs. [2511101 . 
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the Cauchy theorem gives us for the power a^ t (Q 2 ) (y > 0) of this spacelike coupling the following dispersion relation: 

1 f°° da Im a" J -a - ie) 

«p t (Q 2 ) = - / — — (1) 

and r\ — > +0. Eliminating the aforementioned cut sector, and keeping the rest of the discontinuity function Pi (a) 
Im a p t(Q 2 = —a — ie) unchanged, results in the MA/ APT model [TlT - fM] . i.e., the spacelike coupling 

In evaluation of general spacelike observables in pQCD, terms of the type ap t (Q 2 ) In a pt (Q 2 ) \n (Q 2 / 'm 2 ) may appear, 
with v > and (nonnegative) integers k and I. In MA/APT, in principle, the factor ln f (Q 2 /m 2 ) may be included in 
the MA-type analytization (see, for example, Ref. [27]). However, the factor In (Q 2 /m 2 ) is analytic function in the 
Q 2 -complex plane with the exclusion of the nonpositive axis Q 2 < 0. On the other hand, this is also the regime of 
analyticity of Af AA \Q 2 ) [note that Q 2 — is a point of nonanalyticity of A^ AA \q 2 ) because dA^ A \Q 2 ) / dQ 2 = 
oo there]. Therefore, we will not include the factor h/(Q 2 /m 2 ) in the (MA-type) of analytization. With these 
considerations, all the other pQCD terms in such observables, ap t (Q 2 ) hi a p t(Q 2 ), are made analytic in MA/APT by 
the same procedure as a^ t {Q 2 ) in Eq. ([2]) 

(a; t (Q 2 )\n k a pt (Q 2 )) iMA) ^aZ A \Q 2 ) 

\ /an 



i 



Im 

dcr 



a pt (— a — ie) ln fc a p t(— o — ie) 



T Ja=0 (a" + Q 2 ) 



(3) 

In this context, the authors of Refs. [2T | I37K39] derived explicit simple formulas for these MA expressions at one-loop 
level, and the related more involved expressions at higher-loop levels - expressions from which the analyticity structure 
in Q 2 and in v can be more clearly seen. For our calculational purposes, though, the dispersion formulas |3| will be 
applied, numerically, when the MA/ APT model is used. 

In general analytic QCD (anQCD) models, for any spacelike observable T>(Q 2 ) which in pQCD is expanded in 
powers of a pt (p 2 )hi a p t(n 2 ) (where p? ~ Q 2 is the renormalization scale), the analytization method consists in 
principle in a replacement 

a; t (p 2 )\n k a pt (p 2 )^ ( a; t (p 2 ) ln k a pt (p 2 )) = A v , k (p?) . (4) 

V /an 

The explicit procedure of such analytization, in a general anQCD model, is described in Ref. [41] . 

MA/ APT has a relatively specific and interesting property, namely, that the couplings differ from the corresponding 
pQCD couplings in a nonnegligible way even at high \Q 2 \ 

4 MA) (Q 2 ) - a pt (Q 2 ) ~ A 2 /Q 2 (|Q 2 | > A 2 ) . (5) 

As a consequence, a2 AA \Q 2 ) and a p t(Q 2 ) n differ from each other in a nonnegligible way even at high \Q 2 \ [14] ; 
we cannot approximate Af AA \Q 2 ) well with a pt (Q 2 ) even at as high \Q 2 \ as ~ M§. The (only) free parameter in 
MA/APT is the MS scale A, which appears in the usual perturbation expansion of the underlying pQCD coupling 
a p t(Q 2 ; MS) in inverse powers of ln(Q 2 /A 2 ) [i.e., the a p t whose discontinuity functions appear in the relations |2j)-([3j)]. 

It turns out that for MA/APT this scale varies in the interval Ajv /= 5 = 0.260 ± 0.030 GeV. The analysis of high 
energy QCD observables (with \Q\ ~ 10 1 GeV) in Ref. [14] suggests a value of Ajv /= 5 ~ 0.290 GeV [corresponding 
to a s (Mf,MS) « 0.124]. The works of Refs. [131351133 use the value A Nf=5 « 0.260 GeV instead [corresponding 
to a s (M|,MS) w 0.122 and giving the timelike (Minkowskian) coupling 2l 2 (M|) w 0.120]. On the other hand, if 
the "average" e + e~ annihilation ratio R(s) value at y/s — Mz is taken to be « 1.03904 as used, e.g., in Ref. [42], 
MA/APT approach of evaluating R{M 2 ) gives A Nf=5 w 0.225 GeV g3J. 

In comparison, pQCD analyses give for a s (M|) the world average 0.1184 ± 0.0007 [44], corresponding to Aat /=5 = 
0.213 ± 0.008 GeV, i.e., significantly lower scales than for MA/APT. 4 



4 Also the anQCD models of Refs. 15 16 differ nonnegligibly from pQCD, even at high |Q 2 |. 
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Due to the nonnegligible difference, Eq. ([5]), the threshold effects of heavy quarks have to be implemented in 
MA/APT. This is performed by requiring the continuity of the perturbative coupling [whose discontinuities are used 
in MA/APT] 

a pt (Q 2 ) = f Nf (Q 2 /A 2 Nf ), (6) 

at positive "threshold" values Q 2 = Q% lh _^. N > chosen to be for Nf = 4,5,6 the squares of the current heavy 
quark masses m c = = TO5 and mt = me, respectively |14) 

/iV-i^/A^.!) = f N (m% /A%) (N = 4, 5, 6) . (7) 

Further, it is assumed that, for complex Q 2 , the values of a p t(Q 2 ) involve the scale Ajy f determined by the absolute 
value \Q 2 \ 

a P t(Q 2 ) = f N (Q 2 /A 2 N ) for : m 2 N < \Q 2 \ < m 2 N+1 . (8) 

Thus constructed a p t(Q 2 ) gives the piecewise dicontinuous global MA/APT discontinuity functions. For example, for 
Pv'ki' 7 ) = I m [ a pt mfe a pt( — a ~ * e )] 5 this globalization results in 



P^k( a ) = Im(a^(-cr-ie)ln fe a pt (-cr-ie)) 



= Im 



f N (zy\n k f N (z) 



for : (m 2 N < a < m 2 N+1 ) . (9) 

= - CT / A ]V _ie 



Despite the discontinuity of p^(u), the MA/ APT analytic analogs (jij) are analytic functions. 
For the underlying a p t(Q 2 ) = a s (Q 2 )/ir we will use the effective three-loop coupling of the form 

a pt (Q 2 ) = -- r 2 \ , (10) 

where Cj — (3j/Po (j — 1, 2), Q 2 = \Q 2 \ exp(icf>), and W-\ and W+i are the branches of the Lambert function for the 
case < (j> < +7r and — n < (f> < 0, respectively, and 

z(Q 2 ) = ^ exp (-iM/ci) . (11) 



A 2 



We call the scale A 2 (~ 0.1 GeV 2 ) appearing here the Lambert scale. The aforementioned values Aat /= 5 = 0.260±0.030 



GeV correspond to Ajv /= 5 = 0.438 ±0.051 GeV. The solution (10) is a solution to the renormalization group equation 
(RGE) where the beta function has the Pade form f3(a pt ) oc a p r Jl/l](a pt ) 

da P t{Q 2 ) _ _ R 2 [1 + (gi - (c 2 /ci))a pt ] , . 

dlnQ 2 " Poa * [l-(c 2 / Cl )a pt ] ' [U) 

The expansion of beta function /3(a pt ) on the right-hand side gives 

P(a P t) = -Ajflpt (! + c i a pt + c 2 a 2 t +c 3 Op t + . . .) , (13) 

where the higher renormalization scheme parameters Cj (j > 3) are fixed by the leading scheme parameter c 2 : 



c 



C2 /cj . The three-loop scheme parameter c 2 will be fixed so that it is equal to the MS value c 2 (for a given 
Nf). On the other hand, /3 = (1/4)(11 - 2iV>/3) and ci = (3i//3 = (1/4)(102 - 38JV//3)/(ll - 2AT//3) are universal 
constants. The four-loop scheme parameter C3 = PzlPo ls equal here to c 2 /cf, as can be seen from Eq. (13), i.e., the 



value of C3 cannot be adjusted (e.g., to the MS value C3) once the value of c 2 has been fixed. The specific choice 
( 10 1 is convenient because it represents an explicit function 5 and it is thus easy to evaluate it and the corresponding 



5 In Mathematica |48| the Lambert function W n (z) is implemented under the name ProductLog[n, z] 
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discontinuity functions p^(<r), Eq. (jsj), in practice. On the other hand, it approximates the usual (MS) a pt (Q 2 ) to 
three-loops. It has Landau singularities (this is true even for any values of C2). 

Using the aforementioned value A5 = 0.438 ± 0.051 GeV (A5 = 0.260 ± 0.030 GeV), the continuity conditions 
Q give us at other values of N f : A 4 = 0.745 ± 0.077 GeV (A 4 = 0.351 ± 0.036 GeV), A 3 = 0.948 ± 0.085 GeV 
(A 3 = 0.382 ± 0.033 GeV), and A 6 = 0.119 ± 0.015 GeV (A 6 = 0.112 ± 0.014 GeV). The first three quark flavors are 
regarded to be massless. 

In this context, we mention that the elimination of the Landau singularities (i.e., analytization) of a p t(Q 2 ) can be 
performed in various ways, not just in the "minimal" way of Eq. |2]). This analytization, in general results in an 
(analytic) running coupling A\{Q 2 ) which differs from the corresponding perturbative one by a power term 

A 1 (Q 2 )=a pt (Q 2 )+o((^\^j , (Q 2 >A 2 ). (14) 

In the case of MA/ APT, the power index is n = 1, cf. Eq. |5]). This index can be increased to n — 3 jT7l [31] , and 
even to n — 5 [32] , while still maintaining analyticity. In Refs. |29l 130) the problem of finding perturbative analytic 
couplings (in specific schemes) was investigated. It turned out that such couplings always gave significantly too low 
value of the (strangeless) semihadronic tau lepton decay ratio r T , unless the scheme was changed in a drastic way 
which made the perturbation series diverge strongly after the first four terms. 



Therefore, the appearance of the power terms as given in Eq. ( 14 ) seems to be a general feature of procedures 
which eliminate the unphysical (Landau) singularities of the coupling in the complex Q 2 -plane. These power terms 
are of ultraviolet origin and thus contravene the philosophy of the ITEP school [53] . according to which the power 
terms in (inclusive) QCD observables appear in the Operator Product Expansion (OPE) and are of infrared origin. 
Furthermore, in MA/APT we have n = 1 and the terms ~ A 2 /Q 2 appear even in the massless QCD, i.e., in the 
case when such terms are not allowed in the usual pQCD+OPE approach for observables related with the vacuum 
expectation values (such as the Adler function). In general, the ITEP interpretation of OPE must be abandoned in 



those anQCD models in which the index n in Eq. (14) is low (e.g., n < 2), and modified/restricted in others (e.g., 
when n = 3,4). 6 

Despite these reservations, we evaluate in this work quarkonium masses with the MA/APT model of anQCD, in 
order to illustrate the effects of elimination of Landau singularities, within such a context, on the behavior of the 
series. MA/APT model is also relatively simple to apply technically. 

We will also apply another anQCD model to the evaluation of quarkonium masses, namely the so called two-delta 
anQCD |32j . in which the index in Eq. ( |14| is n — 5. Evaluations in such models are technically more demanding, 
though, because the evaluation of the analytic analogs of noninteger powers a pt and of the logarithm-type terms 
a" ln fe a is more involved 1411 . 



B. Two-delta analytic QCD model 



It is possible to construct such anQCD models w hich , at high momenta, practically merge with pQCD; namely, 
with the analytic coupling Ai(Q 2 ) which fulfills Eq. (14) with n > 1. One such model, with n = 3, was constructed 
[T7] 7 Another such anQCD model, which we can call one-delta model, was constructed in Ref. [3T], and also 



Rcf. 



satisfies Eq. (14) with n = 3. The idea was simple: the perturbative discontinuity function p < ^' t \a) = Im a p t(Q 



le 



for a > 0, was kept unchanged for a down to a "pQCD-onset-scale" Mq ~ 1 GeV , and the unknown 
low-energy regime < a < Mq was parametrized with one delta function at a scale M 2 (such that: < M 2 < Mq) 8 



p[ ld \a)=7rf 2 A 2 6(a-M 2 ) 



Q(a - M 2 ) 



(pt)z \ 



(15) 



where for a p t(Q 2 ) [the latter defining p^^icr)] the renormalization scheme with C2 = C3 = . . . = was used. The 
(Lambert) scale parameter A was fixed by the value of q s (M|, MS). The dimensionless parameters Sj = M 2 / A 2 



Nonetheless, even in MA/APT, OPE can be applied, and has successfully been applied - see, for example, the inclusion of higher-twist 
terms in the MA/APT vs pQCD analysis of the Bjorken polarized sum rule at low Q 2 , Refs. [49) . We note that the leading-twist (D = 0) 
term, in such a MA/APT+OPE approach, contains implicitly (parts of) the power term contributions ~ (A 2 /Q 2 ) n , n = 1,2, .... 
The higher couplings A n (Q 2 ) = (a pt (Q 2 ) n ) an were not constructed there, though. 

A similar idea was applied in Refs. 1501 1511 directly to spectral functions of the vector current correlators. Other approaches of 
eliminating unphysical singularities directly from specific (spacelike) obsevables, were presented in Refs. I52H54| . 
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TABLE I: Values of the parameters of the considered 2danQCD model. We consider C2 = —4.76 (Mo = 1.25 GeV) as the 
central representative case. The Lambert scale values in the corresponding cases are for the QCD coupling parameter value 
aP^Mf) =0.1184. 



c 2 = /3 a //3o 


so 


si 


f? 


S2 


/I 


A [GeV] 


M 


-4i(0) 


-2.10 


17.09 


12.523 


0.1815 


0.7796 


0.3462 


0.363 


1.50 


0.544 


-4.76 


23.06 


16.837 


0.2713 


0.8077 


0.5409 


0.260 


1.25 


0.776 


-5.73 


25.01 


18.220 


0.3091 


0.7082 


0.6312 


0.231 


1.15 


1.00 



(j = 0, 1) and f 2 were then fixed by the requirement of n = 3 in Eq. (JT4J) (these are in fact two conditions) and by 
the requirement that the central experimental value of the (V + A) r-decay ratio r T , namely (r T ) exp = 0.203 ± 0.004 
[S51[5S], 9 be reproduced in the model. 

As argued at the end of the previous Subsection, the ITEP interpretation of the OPE requires that index n in 
Eq. (14) be relatively high, 10 4 < n < 9. If, instead of one delta, we use two deltas to parametrize the otherwise 
unknown low-er regime behavior of the spectral function /9i(cr), we obtain an analytic QCD model with Ai(Q 2 ) 
fulfilling the relation Hlh with n = 5, Ref. [35], 



Pi (o-;c 2 ) 



2 

tt Y, /JA 2 6(<r - M?) + B(a - Ml) x p^\a; c 2 ) 

2 



)+e( S - S „)x r « pt) ( S ;c 2 ) 



(16) 



(17) 



Here, the dimensionless squared energy s = a /A 2 was introduced, and five dimensionless parameters: sj = M?/A 

(j = 0,1,2) and f 2 (k = 1,2). Further, c 2 = is the scheme parameter, and r^ pt '(s;c 2 ) = p < f t \a;c 2 ) = 

Im a p t(Q 2 = —o — ie; c 2 ) is the perturbativ e sp e ctra l function of a P t (Q 2 ) in t erms of the dimensionless energy variable 



s = a/A , where a p t(Q ) is given by Eqs. (10l-(ll) [see also Eqs. ( 12 )-( 13 )] . The spacelike coupling A\(Q ) is then 
obtained by the dispersion relation 



4 2d) (Q 2 ) 



(2d)/ \ 

Pi W) 



./.T-o da ( CT + Q 2 ) ^[( u + s j) ' 77 



E 



f 2 



„(pt) 



0;c 2 ) 



ds - 



(18) 



where we denoted u = Q 2 /A 2 . This two-delta analytic QCD model (2danQCD) is to be applied in the low-momentum 
where n/ = 3 (considering the current masses of u, d and s zero): |Q 2 | < (2m c ) 2 (« 6.45 GeV 2 ). At 



> (2m c ) (where rif > 4) it is replaced by the underlying pQCD of the coupling a pt (Q ), Eq. (10 1 



regime, i.e. 
higher \Q 2 \ 

The Lambert scale A 2 = A 2 lf=3 (~ 0.1 GeV 2 ), which is the only dimensional parameter of the model, is determined 

by the world average value a pt (Af§;MS) = (0.1184 ± 0.0007)/7r, 44J. For example, when c 2 = -4.76, we obtain 
A 1=3 (0.260 ± 0.008) GeV. On the other hand, the five aforementioned dimensionless parameters Sj and f 2 are fixed, 
independently of the value of A, by the condition Eq. ( 14 ) with n = 5 (these are four conditions) and by the condition 
that the aforementioned central experimental value of the (V + A) r-decay ratio r T be reproduced by the model, 



namely r T = 0.203. More specifically, the (four) conditions Eq. (14) (with n = 5) determine the four relations Sj(so) 
and f 2 (so) {j = 1,2), and the value of s (= M 2 /A 2 ) is then determined (if A is already fixed) by the condition 
r T (V + A) — 0.203. The (scheme) parameter c 2 remains a free parameter. Varying c 2 results in variation of the 
pQCD-onset scale Mq (= ^/sqA) and in variation of Ai(Q 2 ) at low \Q 2 \. For phenomenological reasons, we wish to 
have Mo below the r-lepton mass m T 1.78 GeV) and to have the largest |„4i(<3 2 )| < 1.0. More specifically, we 
vary c 2 in such a way that 1.0 GeV < M < 1.5 GeV and -4i(0) < 1.0. This results in the central value c 2 = —4.46 
for which M = 1.25 GeV (and ^i(0) = 0.776), and two border values: c x = -2.10 for which M = 1.50 GeV (and 
^(0) = 0.544); c 2 = -5.73 for which A\{$) = 1.00 (and M = 1.15 GeV). 11 The model parameters for the central 



9 r T is the QCD part of the V+A decay ratio R T (AS = 0), with the (small) quark mass effects subtracted. It is normalized in the 
canonical way: r T = a + 0(a 2 ). 

Small-size instanton effects may lead to the conditions Eq. | |l4| to be valid only for up to index n such that In is the largest dimension 
of condensates not affected by the instantons. Instanton-antiinstanton gas scenarios lead to n < 4/?o (= 9 for ne = 3), cf. Ref. 1351 . 

1 In Ref. |32| we included the case C2 = —7.15, for which Mo = 1.0 GeV. However, in that case, .4i(0) = 2.29, and this indicates that 
the model has infrared instabilities at such low C2 values. 
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and the two border cases are given in Table [Tj 

Variation of the scheme parameter c%, at fixed central world average value q p t(M§; MS) = 0.1184/7T, implied the 
corresponding variation of the parameters Sj, /| and A in Table |lj such that the relation (14) with n — 5 is fulfilled 
and r-decay ratio value r T = 0.203 is reproduced. 

On the other hand, if we vary the QCD coupling parameter within the world average interval, a p t(M§; MS) = 
(0.1184 ± 0.0007)/7r, [44], while keeping C2 (= —4.76) and all the dimensionless parameters fixed, then only the scale 



parameter A varies, A = 0.260± 0.008 GeV, reflecting the changing a pt (M§; MS). The fulfillment of the relation (1 14h 
with n = 5, at given C2, fixes the values of the four dimensionless parameters Sk and /? (k = 1,2) in terms of the 
parameter so (= 23.06). The resulting value of r-decay ratio then varies with a p t(M§; MS): r T — 0.203 ± 0.006. This 
is acceptably compatible with the measured values (r T ) cxp = 0.203 ± 0.004 55, ?55] . 
For further details on 2danQCD model we refer to Ref. |32j . 

The higher order couplings A v (Q 2 ) (analogs of a^ t (Q 2 ) ) in this model are then obtained according to the construction 
for general anQCD models, explained in Refs. [20l [21] for integer v, and Ref. [41] for general (noninteger) v. Below 
we will summarize the basic aspects of such construction. 

Since the general anQCD models, with the exception of MA/ APT, have at low u (< 1 GeV 2 ) different discontinuity 
function than the pQCD coupling a pt (Q 2 ), we cannot use the MA/APT method, Eq. for the construction of 
the analytic analogs of a pt ln fe a pt - The analogs of the integer powers a™ t in such general models were constructed in 
Refs. [201 [H], where it was shown that it is imperative to construct first the analogs of the logarithmic derivatives of 
a p t in the following way: 12 



fd k a pt (Q 2 )\ _ d k A 1 {Q 2 ) 



\8(\nQ 2 ) k J an d(\nQ 
In pQCD, the logarithmic derivatives 



(* = 0,1,2,...). (19) 



<* = °. (20) 

are related with the powers of a p t = a s /f in the following way [using renormalization group equations (RGE) in 
pQCD]: 

a P t,2 = a p t + cia pt + c 2a pt + c 3« P t ^ » ( 21 ) 



Opt,3 = < + ^c iapt + Uc 2 + ^c 2 \ a pt + • • • , (22) 



13 

a P t,4 = a pt + y ciapt + ■ • ■ , a pt , 5 = apt H , etc. (23) 



This means that the powers of a pt are linear combinations of logarithmic derivatives 

2 ~ ~ /5 2 \ . / 28 , 22 \ . 
a t = a pt ,2 - cia pt ,3 + x c i ~ c 2 a pt ,4 + - — q + — cic 2 - c 3 a pt ,5 H , (24) 



2 , .j ^ 3 x 3 

a P t = « P t,3 - ^cia p t,4 + (^y c i ~ 3c 2^) «pt,5 H , (25) 

13 

a pt = a pt,4 ^-cia pt: 5 H , a p t = a P t,5 H , etc. (26) 



These relations, in conjunction with the analytization Eq. (19), imply that the analytic analogs Ak of powers a k t , in 
general anQCD models, can be expressed as linear combinations of the logarithmic derivatives 

Ak+lip ] (3 k k\ d(^f ' (fc- 1 ' 2 ---)- (27) 



12 If the analytization is performed in any other way, the renormalization scale and scheme dependence of the resulting truncated analytic 
series of any observable T>(Q 2 ) will in general increase (instead of decrease) when the number of terms in the series increases, cf. |20II21I . 



in the following form: 



A 2 


= A 2 


A 3 


= A 3 


A 4 


= A 4 



c%A 3 



c 2 Aa 



2N 



22 



5 
2 

13 

y 



c\Ai 4 
ciA 5 



28 



3c 2 A 5 



A 5 = A 5 



-cic 2 - c 3 A 



etc. 



In Ref. 



this analytization was extended to the case when k ^ v is noninteger 13 

1 (-1) f°° da 



A u+ i(Q z 



7T + 1) J 



-pi{<j)lA_ v 



Q' 2 



(-1 < v) , 



(28) 
(29) 
(30) 

(31) 



where, as always, pi(<r) = Im A\(Q 2 = —a — ie), and Li_„(z) is the poly logarithm function. 14 The corresponding 
analogs of powers a pt are then obtained by using the general relations 

oo 



pt 



(32) 



and the linearity of analytization, i.e. 



AAQ 2 



*p t (Q 2 )) an = A V {Q 2 ) + ]T k m (y) A„ +m {Q 2 ) iy > -i) 



(33) 



m>l 



The cocfhcicnts k m {v), for general real v and positive integer m, were calculated in Ref. [41] ■ and are combinations 
of Gamma functions and its derivatives, with arguments involving v + i [I being various integers). 
Furthermore, since a v ln fc a = d k a u jdv k , the linearity of analytization then implies 



A v . k {Q 2 )^ (^ t (Q 2 )ln fe a pt (Q 2 )) i 



d k A„(Q 2 ) 
dv k 

d k A v (Q 2 ) 
dv k 



EE 

m>l s=0 



d s k m {v) d k ~ s A u+m (Q 2 ) 



dv k ~ 



(y > -1) , (34) 



where in the terms on the right-hand side we use expressions for A^+ m (Q 2 ) obtained by Eq. (31 1 . Comparing Eqs. (34) 



and (31 ), we see that the terms in the above sum represent integrals over the scale a involving the basic discontinuity 
function of the model (pi) and derivatives of the polylogarithm function Li_j,_ m+1 with respect to its index v. In the 
evaluation of the binding energy E q q we will encounter the logarithmic terms of the type (34) with v integer [y = n); 



however, the derivatives with respect to index v in Eq. (34) imply that we must know the behavior of A v around the 



integer value v — n, i.e., we need to use here the expressions ( |31[ ) for noninteger or a version of it with improved 
integration convergence, Eq. (22) of Ref. [H] (cf. also the earlier footnote 13). 



The analytization of expansions involving powers a^ t (p 2 ) and ap t (/u 2 ) In a pt (^ 2 ) now consists simply in the replace- 



a^ 2 ) ^ (a» t {p 2 )) an =A„{p 2 ) , 
a; t (n 2 )ln k a pt (^ 2 ) ^ U t {p 2 ) ln fc a pt (p 2 )) = A^) 

\ /an 



ments 



(35) 

In the special case of MA/APT, this analytization is, in principle, equivalent to that of Eqs. ([2|-([3|, as argued in 
Refs. [HJUT]. The (small) differences between the two analytizations in MA/APT arise due to somewhat different 
truncations applied in the two approaches. One is the truncation in the loop-expansion (< ap^ 1 ) usually applied 
to the running coupling a p t(/i 2 ), which is then reflected in the spectral functions Ima pt (— a — ie) and Im[ap t (— a — 



13 This relation was also reformulated so as to be applicable in a larger v interval: —2 < v, cf.Eq. (22) of Ref. | 4 1 1 . 

14 In Mathematica |48| it is implemented under the name PolyLog[— u, z]. 
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FIG. 1: The analytic couplings Ai(fJ- 2 ) at low positive fi 2 (= Q 2 ), for MA/APT and 2danQCD. Included are also pQCD couplings a pt (fj, 2 ) 
in MS and Lambert scheme. See the text for details. 
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\ \ « p , 2 (MS) 




" A 2 t2d \ \ 
















A (MA) ^^^o . 
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FIG. 2: (a) The squares of the pQCD couplings, a pt , and their analytic analogs Ai, at low positive fi ; (b) the logarithmic 
derivatives a p t,2 and A2 at low positive fi 2 . See the text for details. 



ie) ln fc a pt (— a — ie)] in Eqs. ([^-([sj) • ln the general approach, Eqs. (31 1 - ( 34 ) , which can be applied also to the MA/APT 



a) = Im a pt (Q = —& — ie) in the 



analytic model, the same kind of truncation for a pt (^ 2 ) can be applied in p^f^l 

integral (I31J), and the corresponding loop-truncation m < m max = I — v + 1 in the summations (33) and (34). 

]we present the analytic running couplings _4i(/i 2 ) of 2danQCD and MA/APT, at positive low values of ji 2 

-4.471, 

^■j — u 2 / °i ~ f° r 3 — 3; in 2danQCD: see the second line of Table Q <■■ 1.76, r, '-'■•■".'-^ 



In Fig. 

(= Q 2 ), for the central input values. In MA/APT this means: A„ /=3 = 0.948 GeV, c 2 = c 2 (MS, Nf = 3) 
Cj = 4 /c[~ 2 for j > 3; in 2danQCD: see the second line of Table c 2 = -4.76, Cj = qfVcr 2 for 

j > 3. For 

comparison, we also include the pQCD couplings a pt (fi 2 ), in the Lambert renormalization scheme (i.e., the scheme of 
2danQCD), and in the 4-loop MS scheme: c 2 = c 2 (MS, N f = 3) +4.471; c 3 = c 3 (MS, iV/ = 3) « +20.99; 9 = for 
j > 4. In Fig. [2j(a) we present the higher couplings >4 2 and a pt , where in MA/APT the coupling A2 was constructed 

according to Eq. ([2]), and in 2danQCD according to Eq. (28) truncated at A5 term [cf. also Eq. (27)]. While we can 
see in Fig. [I] that the Lambert a p t and 2danQCD A\ practically agree with each other down to scales /1 2 ~ 1 GeV 2 , 
we can also see in Fig. [2^a) that Lambert a 2 t differs more significantly from the constructed 2danQCD analog A2 
even at somewhat higher scales, this difference being primarily the consequence of the aforementioned truncation 
at ^5 in the construction of A2. In Fig. [^b) we present, for comparison, the logarithmic derivatives a p t. 2 (/i 2 ) and 
^(/J?) [cf. Eqs. 20) and (27)], where we can see that the difference between these couplings in 2danQCD and the 



corresponding (Lambert) pQCD is again negligible for fj, 2 > 1 GeV^. 

The 2danQCD model has been successfully applied [33J in the analysis of (V - 
with OPE. 



■ A) r-decay data via Borel sum rules 



10 



III. PERTURBATION EXPANSION FOR HEAVY qq GROUND STATE ENERGY 



The analysis of nonrelativistic potential is the starting point for the determination of the ground state energy of 
qq and thus of the mass of such systems. The main input in these calculations is the mass m g (m 2 ) (here denoted 
simply as m q ), also called MS quark mass. The masses of the heavy quarkonia qq (T(IS') when q = b, J/ip(lS) when 
q = c) are well measured, and this allows us to extract the corresponding mass fn q . By evaluating an observable, 
such as the quark-antiquark binding energy here, within anQCD models, at least part of the (chirality-conserving) 
nonperturbative effects gets included in the leading-twist term via the analytization, such as Eqs. ([2|-([3]) in MA/APT, 
and (35 1 in general analytic QCD models. 

The coefficients in the (leading- twist) perturbation expansion of the ground state binding energy E qq (i.e., with: 
n = 1 and £ = Q) of heavy quarkonium qq in powers of et pt were obtained up to all terms 0(a pt ) in Ref. [58] . the 
terms 0(a pt lna pt ) in Ref. [59], and all terms up to 0(a pt ) (including logarithmic) are given in Ref. [7]. The last term 
(~ Op t ) is now completely known since the parameter 03 from the static potential is now known, Refs. [5J[B]. The 
general structure of the (leading-twist term of the) ground state binding energy in pQCD is 



p(pt) 



"^^ 2 «pt(M 2 ){l + a P t(M 2 ) [#1,0 + K hl L pt {fi 2 )] + a 2 pt {fi 2 ) [K 2>0 + K^L^p 2 ) + K 2t2 L 2 pt {^ 2 )] 
-a^p 2 ) [#3,0,0 + #3,0,1 hia pt (A* 2 ) + K 3A L pt (^ 2 ) + K 3 , 2 L 2 pt (^ 2 ) + K 3i3 L 3 pt (^ 2 )] + 0(a%)} , (36) 



where 



>p1 



(M ;rn 



In 



((47r/3)m 9 )2 a 2 t(/i 2 



(37) 



/i 2 is the (square of the) renormalization scale, m q is the pole mass of the quark, and the coefficients A"j,fc can be 
obtained by combining the results of the aforementioned literature, presented also in Appendix [X] The typical scale 
of the process is a so called soft scale Q 2 (= — q 2 ) which represents a typical quark-antiquark momentum transfer 
inside the quarkonium (Q 2 ~ m 2 a 2 ). The "soft" renormalization scale /x 2 = // 2 is expected to be of the same order 



(/s~l;Q2~m?a2) 



(38) 



The quarkonium mass is then 



M qq = 2m q + E qq {m q ) 



(39) 



In principle, the input quantity here could be the quark pole mass m q . However, this mass m q , in contrast to the 
mass fn q , suffers from the strong infrared renormalon ambiguity (at Borel parameter value b— 1/2, =>- 5m q ~ Aqcd)- 



This ambiguity must cancel in the physical sum, Eq. ( 39 ) 



Therefore, it is more convenient to use as input a renormalon ambiguity- free input mass, such as m q , and we will 
follow this approach. 

In the case of the bottom quark, before we relate the pole mass nib with the mass m&, at renormalization energies 
(i = WL^ where Nf = 4 (and where the charm quark mass is considered zero, i.e., decoupled), the effects Smi, of the 
nonzero massa m c ^ have to be subtracted, and they are |60j 



(5m b ) mc = (5m b )£] + (6m b )<£ « 0.025 ± 0.005 GeV 



(40) 



These effects were calculated in Ref. [SD] in pQCD, at the hard renormalization scale /1 2 = m 2 . We checked that they 
do not get significantly modified in MA/APT and in 2danQCD. 

The pole mass m q and the mass m q are then related via the relation 



TUX — $7Tl 

1 = l + R (a pt (m 2 ) + r ia 2 t (m 2 ) + r 2 a 3 pt (m 2 ) + r 3 a 4 pt (m 2 q )) + 0{a b pt ) , 



(41) 



where Sm q is zero when q = c, and is given by Eq. ( 40 ) when q = b 



5 mq ^\^ mb) ^ {q = b) 



(Q = c) 



(42) 
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The coefficients Ro, ft, and r 2 are known (Refs. [51] , [52] . [6"3"rl6"5] . respectively) 

i?o = ^ , i? ri = 6.248/3 - 3.739 , 

R r 2 = 23.497/^ + 6.248/3i + 1.019/3 - 29.94 , 



(43) 
(44) 



where /3q = (11 — 2Nf/3)/4 and /3\ = (102 — 38JV//3)/16 are the two universal aforementioned RGE coefficients, and 
Nf is the number of flavors of quarks lighter than q. Specifically, the values are: r\ — 7.739 and r 2 — 87.224 for 
Nf = 3; and r\ = 6.958 and r 2 = 70.659 for Nf — 4. On the other hand, in Appendix [B] we estimate the values of r 3 
to a reasonably high level of precision (with less than 4 % uncertainty) by a method which uses the structure of the 
leading infrared renormalon (at b = 1/2) of the quantity m q /m q 



r 3 (N f = 3) = 1339.4 , r 3 (N f = 4) = 987.3 



(45) 



While the relation ( 41 ) is written here at the "hard" renormalization scale /x 2 = m 2 , it is straightforward to reexpress 



the sum on the right-hand side of Eq. (41 ) at a different, lower, scale \i . In Appendix \B\ this reexpression is presented 



explicitly, under the assumption that during the lowering of the scale, fn q — > (j, , we do not cross quark threshold. The 
goal is to express in the perturbation expansion of the binding energy (36), where the renormalization scale is soft, 
the pole mass m q via the mass m q , and for this we need the relation (41 1 at the soft renormalization scale. It turns 
out that for the bb system the hard scale is mj ss 4 GeV, i.e., the scale where Nf = 4; and the soft scale, Eq. (38), is 
fj, 3 w 2 GeV, i.e., the scale where it is more reasonable to expect Nf = 3. 15 In this case, we take into account also 
the (three-loop) quark threshold transition Nf = 4 H> 3 at fi 2 — (2m c ) 2 , Ref. [66]. We thus obtain the relation (41) 
reexpressed at the soft scale /z 2 of Eq. ( 38 1 



m q (ml; fi 2 ) = 8m q + m q I 1 + -a pt (^) [l + a pt (^ s ) ri {fi 2 s ) + al t {fi 2 s )r 2 (p 2 s ) + a^ t (/x 2 )r 3 (/z 2 )] 



(46) 



Further, the renormalization scheme can also be varied in this relation, i.e., the changes of the scheme parameters 
Cj = /3j//3o (j = 2,3) from the usual MS scheme to other schemes affect correspondingly the values of the coupling 
a p t(/i 2 ) and of the coefficients rfc(/u 2 ). The relationship between a p t's at two different scales and in two different 
renormalization schemes is summarized in Appendix [Cj where we also summarize the (3-loop) connection of a p t 's 
across the quark threshold. 

After carefully performing all these transformations, we can rewrite the original expansion ( |36[ ) for the binding 
energy in terms of the m q quark mass, at any soft renormalization scale /x s and in any chosen renormalization scheme 
(c 2 ,c 3 , • . •) 



\Q 2 ;m 2 q ;N f =3) = -g{m q + 5m 3 )7r 2 a^ t (/i 2 )|l + a pt (p 2 ) [JC l<0 + /Ci,i£ pt (/x 2 )] 
+Op t (Ai 2 ) [£2,0 + £2,i£ P t(^ 2 ) + /C 2 ,2£pt(M 2 )] 

+a pt (/x 2 ) [£3,0,0 + £3,0,1 lna pt (/i 2 ) + £3,1^^) + ^C 2 ^) + /C 3 , 3 /: pt (/i 2 )] + 0«)} , (47) 
where the logarithm contains now fn q quark mass 



£ P t(M 2 ;m 2 ) = o ln 



[(47r/3)(m q + Sm q )] 2 a 2 pt ^ 2 ) / 



(48) 



mass m q , but rather 



We note that the new (renormalon-ambiguity-free) mass which appears naturally in this expansion is not exactly the 

(49) 



• ' Q 1 fn c (q = c) 



where (5mi,) mc is given by Eq. (40) when q = b. 



15 Usually, the quark thresholds are taken at 2m,. In the case of Nf = 4 1 — > 3 transition, this is about 2.5 GeV, above the soft scale of 
the bb system. 
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The aforementioned soft "process scale" Q 2 (~ m q a 2 ) can be regarded, at least formally, to be a variable complex 
scale. Therefore, the binding energy E qq (Q 2 ;m 2 ) is, formally, a spacelike observable analytic in Q 2 . The analytization 
of E„f\Q 2 s -ml) of Eq. |47l), according to Eqs. S [or, in MA/APT: Eqs. (pl-Jil], then leads to 



E q f(Q 2 ;m 2 q ) = - -m q 7r 2 {a 2 (kQ 2 s ) + [JC 1<0 A 3 (kQ 2 s ) + £1,183,1 («<#)] 
+ [JC 2fi Ai{nQ 2 s ) + JC 2A B^{kQ 2 s ) + £2,2-64,2 faQ*)] 

+ [/Ca.o.oA (kQ 2 s ) + £3,0,1^5,1 (kQD + /C 3) iB Bl i(«Q2) + £3,2S 5 , 2 (kQ') + A :i .;5'r, .ji/.-Q;'; + G(A, 4 )} , (50) 

where we use, for simplicity, the notation of Eq. (33) with Eq. (27) for A„'s (y = k integer now) [in MA/ APT: Eq. Q], 
and denote by B n + 2 ,j the following: 



ft 



n+2,j 



where 



s=Q 



2J 



((4n/3)m q ) 2 a 2 pt (KQ 2 ) 



(-2)* f- s (nQ 2 s )An+2A^Ql) 



(51) 



/=ln 



(4ir/3) 2 M 2 



(52) 



and A n+2 ,s were defined in Eq. (34 1 [in MA/APT: Eq. }3l]. 



Here, we have j 
freedom. 



1, . . . ,n, implying that B n + 2 j(KQ2) ~^ faster than <Zp t (K(5s) when \Q 2 \ — > 00, by asymptotic 



The relation ( 46 ) in its analytized form is 

~m q (ml; n 2 s ) = 5m q 



fe g + m j 1 + ^ [ii(^) + ri(M^2(M») + r 2 (^)i3(^) + r 3 (M^^)] +0(A 5 )} . (53) 

Since this relation includes at its highest order the term ~ A4 (~ a£t), m general analytic QCD it is consistent to 
use in Eq. ( |53"| for the expressions A n (n — 1,2,3,4) those given in Eq. (33) [or: Eqs. (28)-(30)] with the sum there 
truncated at (and including) A4 (~ A4 ~ Opt)- On the other hand, the expression (501 for the binding energy includes 
at its highest order the terms ~ A5 (more precisely: ~ ^5,3), therefore it is consistent to use there for the expressions 
A n (2 < n < 5) those of Eq. (33) [or: Eqs. ( 28 )-( 30 )] with the sum truncated at (and including) A5, and expressions 
A n ,k of Eq. (34) truncated at (and including) d k A v +^/dv k \ u=0 (~ a™ t ln fc a pt ), where k = 1, 2, 3. 



The aforementioned soft scale Q 2 (~ m 2 q a 2 ) can be fixed in pQCD, by convention, by the condition that all the 
logarithmic terms in the expansion ( 47 ) disappear 



Qs,pt 



(4ir/3) 2 m 2 q al t (Q s ^j 



(54) 



where fh q is defined in Eq. (49). This type of condition, in analytic QCD, would correspond to fixing a referential soft 
scale Q s by requiring B n +ij(Q 2 ) = 0, for various n's and j = 1, . . . , n, meaning that the fixing is not unique since 
it depends on n and j. Our convention will be to require that the leading logarithmic term i3„ +lj in Eq. ( |50| , i.e., 
63,1, be zero at such scale 



B 3 ,i(Q2) = 



(55) 



It will turn out that this condition has a solution in the case of bb (T(1S)) in 2danQCD, but not in J/ip(lS) in that 
model, and not in any case of the MA/APT model. In such respective cases, we will use simply the following, more 
naive, analogs of the pQCD condition |54| ): 

(56) 
(57) 



Q 2 = (Ai:/3) 2 m 2 q A 2 (Q 2 s ) (J/ip{lS) in 2danQCD) 
(MA) = (4^/3) 2 m^4 MA) (Q^ (MA) ) (MA/APT) . 



These measures of the typical momentum scale of the (nonrelativistic) quark inside the quarkonium are rather low, 
sa 2 GeV in T(15), and < 1 GeV in J/ijj(lS). In pQCD such scales are problematic, because they are not far 
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away from the unphysical (Landau) singularities of a p t(Q 2 ); in ananlytic QCD models, no such problems appear in 
principle. 

We recall that the pole mass m q and the static potential V(r) both contain the leading infrared renormalon (b = 1/2) 
singularities, and that cancellation of these singularities takes place in the sum 2m q + V(r) As a consequence, 

this cancellation must take place also in the quarkonium mass 2m q + E qq [60, 70 74 , more specifically, in the sum 
2m q + E qq (s) where E qq (s) = (1| V(t-) |1) (~ V(r s )) is the soft part of the binding energy E qq , and |1) denotes the 
(ground) state of the quarkonium. The typical soft distance r s in the quarkonium is ~ l/-v/Qf where Q 2 is the 
soft scale of Eqs. ( 55 1-( 57 ) , i.e. l/r s ~ (m q ira) (where a — a p t or A\). However, since V(r s ) oc l/r s , we have 
E qq (s)/m q ~ V (r s) / m q oc l/(r s m q ) ~ (7ra), this thus leading to the so called "power mismatch" in the renormalon 
cancellation in the pQCD expansions 2m q + E qq (s) [75] (see also: [70]): the terms ~ a" in 2m q /m q tend to cancel 



pt 



it is 



numerically the terms ~ o" +1 in E qq /m q . Therefore, since the binding energy E qq is now known up to 
very convenient to have the relation m q /m q up to ~ a 4 , i.e., to have a good estimate for the coefficient r% and to use 
it, so that the effects of renormalon cancellation in 2m q + E qq (s) can be seen numerically more clearly. Of course, 
this cancellation, term-by-term, should be numerically more precise (at sufficiently high orders) if the renormalization 



scales fi s used in m q /m q and in E qq (s)/m q , Eqs. (461 and (47 1, are taken to be equal. This renormalon cancellation 



will be our guiding principle for the separation of the ultrasoft (us) part from the binding energy 



E qq 



E qq (s) + E qq (us) . 
■ 2 and the us part of the binding energy is 



(58) 



Typical us scales are fi us ~ m q a s , and the us part of the binding energy is ~ m q a b In a. We can parametrize the s-us 
separation by a dimensionless parameter k 8 / us such that the s-us factorization scale fj,f is written as 



lqCts(Ql 



\3/2 



/us 



(Qs 



a/2 



(59) 



It is expected that usually k s / us ~ 1, but it does not have to be so always. The us part can be rewritten, in terms of 
k s /us as (cf. Ref. [74] ) 



E q f(us) 



m 9 
and in analytic QCD correspondingly 



m q TT 2 [£3,0,0 (us)al t (nl s ) + JC 3 ,o,i( us ) a lt (mL) 111 a P t (mL) + °(«pt)] 



(60) 



E q f(us) 



-m q n 2 [K,3,ofi(us)A 5 {^l s ) + /C 3i0 ,i( us )A,i(A i L) + , 



(61) 



where fi us is a us renormalization scale (fi us ~ m q a 2 ), and the two us coefficients are 

£3,o,o(ws) = [27.512 + 7.098 In tt - 14.196 \n(k s/us )] n 3 , 
fC 3 , ,i{us) = 7.098tt 3 . 



(62) 
(63) 



The expansion of the soft (s) part E qq (s) of the binding energy is then, according to Eq. (58), the same as the 
expansions (47) and (50), with the exception of the replacements 16 of two coefficients /C3 0,0 and K3 1 



EqT' (s) = E ( q f' an) (/C 3 ,o,o H' /C 3 ,o,o - £3,0, o(ws); £3,0,1 n- £3,0,1 - £3,0,1 (■"«)) 



(64) 



The s-us factorization, i.e., the parameter k s / usi will then be determined, in each model, by requiring that the leading 
infrared renormalon cancellation in 2m q + E qq (s) be exact at the last available order, i.e., that the C(a 4 ) term in 
2m q (Q 2 ,fn 2 ) and the term C(a 5 ) 17 in E qq (s;Ql;fn 2 ) cancel exactly. 

The us part of the quarkonium mass, E qq (us; Q^ s ; rn 2 ), will be evaluated in each case according to a procedure 
which takes into account those problems of low-scale evaluations which appear in the considered model (pQCD, 
MA/APT, 2danQCD). 



The coefficients ^C3,o,o( u *) an d fc3,o,i(us), representing the leading part of the (quasi) observable E qq (us), are renormalization scale 
(fius) and scheme independent. 

The latter term includes all 0(a 5 In* a) terms (k = 0, 1, 2, 3). 
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IV. NUMERICAL RESULTS 

A. Bottom mass extraction 

In this Section, we extract from the mass of bb quarkonium the mass = ra^mj) in MA/APT, 2danQCD model, 
and in pQCD in two renormalization schemes (MS and in the Lambert scheme of 2danQCD). For this, we use the 
relation between mj, and the well-measured mass of the bb quarkonium T(1S) [U] 

2m b (mt;[£)+E bl (s;Ql-,ml[i 2 s ,N f = 3) + E bb (us; Q 2 us ;M 2 b ; tfj = M^(= 9.460 GeV) . (65) 

The dependence on the (soft) renormalization scale a 2 in the pole mass m b and in the soft binding energy E bb (s) 
occurs due to the truncation of the series for these two quantities. For the same reason, the ultrasoft binding energy 
E bb (us) has strong dependence on the ultrasoft renormalization scale fj,^ due to the drastic truncation of this quantity 
at its leading order (~ a 5 In a). As mentioned in the previous Section, the separation of the s and us parts of the 



binding energy will be performed here by determining the s-us separation parameter k s / us , Eqs. (58)- (64 1, by the 
requirement of cancellation of the leading renormalon in 2m b + E bb {s). 

In contrast to the other three models, MA/ APT gives a very small central value for the s-us separation parameter 
ks/us ~ 10. ~ 12 . This reflects the difficulty in the MA/APT scenario to exactly enforce the leading renormalon 
cancellation of the ~ term of 2m b with the corresponding ~ A§ term of E bb (s). If, on the other hand, we impose 
in MA/ APT the condition k s / us ~ 1, more specifically the central value k s / us = 1 and variation of that parameter in 
the interval (0.1, 10.), the results change somewhat, the central extracted value of fn b increases by about 0.050 GeV, 
and the absolute values of the us part of the binding energy and of various other uncertainties of m b get significantly 
reduced. We will consider in MA/APT only the natural range ~ 1 of the s-us separation parameter, k s / us = l.O^o'g, 
rather than the exceedingly small values of k s / us required by the exact renormalon cancellation. In the other three 
models (2danQCD, and in Lambert and MS pQCD), the renormalon cancellation is imposed without any problems, 
resulting in the values of k s / us within the interval between 10 _1 and 10 1 . 

In the previous Section we mentioned that we take Nf = 3 for the number of active flavors in the binding energy, 
i.e., in this case the m c mass is considered to be infinite (decoupled). It turns out that, while the effects of the 



finiteness of m c cannot be neglected in the relation between mb and Eqs. (40l-(41|, these effects can be safely 
neglected in the binding energy E bb {Nf = 3) , cf. Ref. [60 based on Refs. [6"2l 1761 177] , 

Application of the formalism described in Subsec. II A and Subsec. II B for the calculation of the couplings of the 



analytic QCD models MA/APT and 2danQCD, and in Sec. Ill for the calculation of 2m b , E bb (s) and E bb (us) in terms 
of these couplings, then gives us the following results: 

m 6 (MA/APT) = (4.165±0.002 (u / + [Un> ' l ' -(U)J{) 



+ (;o:2o3)^ ) TO - 005( - ) } GeV (66) 

4.165 ± 0.022 GeV , with : Q S , (M A) = 1-60 GeV , k s/us = 1.0 . (67) 



m b (2danQCD) = |4.353 + ( + °;°^ ) f + ( _Joi6 ) ( ; } T °- 005 K) 

T0.005 (mc) | GeV (68) 



{us) V / (s/us) 

-0.023 \ , / +0.017 



-0.034 , ' \ -0.025 . , 

7 (ca) ^ ' (Mi) 

= 4.353 ± 0.084 GeV , with : Q s = 2.08 GeV , k s/us = 0.241 . (69) 

We will comment on the above uncertainties below. For completeness, we give here also the results of the same kind 
of analysis in pQCD, first in the Lambert renormalization scheme (i.e., the scheme as used in 2danQCD: C2 = —4.76, 
Cj = c^ 1 /c[~ 2 for j > 3); and in (4-loop) MS scheme: 

m 6 (pQCDLamb.) = [4.382 + ( ) + ( ~° Q ™ ) ± 0.010 (a#) 

L \ ' / (us) V ' / (sjus) 

iZ) l +(:2Z) r ^°- oo5 (-)} GeV (70) 

/ (c 2 ) V / (Mi) J 

4.382 ± 0.111 GeV , with : Q a<pt = 1.73 GeV , k s/us = 0.306 . (71) 
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k \ / (us) \ / {s/USj \ / (CKsJ 

laW^^^^V^GeV (72) 
= 4.505 ± 0.231 GeV , with : Q s , pt = 1.87 GeV , k s/us = 0.248 . (73) 

The value of the s-us separation parameter k s / us was determined in all cases by the aforementioned renormalon 
cancellation in the sum Mx(is)(s) = 2m;, + E bb (s;k s / us ), except in the MA/APT case, as discussed above. Below 
we present these resulting sums, for the central choices of the aforementioned four results, where we combine in each 
parenthesis the (positive) terms ~ a n of 2m b and the corresponding (negative) terms ~ a n+1 of E bb (s) (n = 0, 1, 2, 3, 4), 
in order to see more clearly the tendency of the renormalon cancellation; the us part is given separately 

M T(15) (s; MA/APT) = 8.380 + (1.151 - 0.152) + (0.281 - 0.173) + (0.137 - 0.152) + (0.075 - 0.092) GeV 

= 8.380 + 0.999 + 0.107 - 0.015 - 0.017 GeV (= 9.454 GeV) , (74) 
E bb (us; MA/APT) = 0.006 =p 0.003 GeV , (m b = 4.165 GeV, k s/us = 1.0) ; (75) 

M T(1S) (s; 2danQCD) = 8.756 + (0.999 - 0.131) + (0.373 - 0.222) + (0.233 - 0.193) + (0.568 - 0.568) GeV 

= 8.756 + 0.868 + 0.151 + 0.040 + 0.000 GeV (= 9.814 GeV) , (76) 
E b - b (us; 2danQCD) = -0.353 ± 0.151 GeV , (m b = 4.353 GeV, k s/us = 0.241) ; (77) 

M r{1S ) (s; pQCDLamb.) = 8.814 + (1.095 - 0.170) + (0.327 - 0.220) + (0.346 - 0.376) + (0.413 - 0.413) GeV 

= 8.814 + 0.925 + 0.107 - 0.031 + 0.000 GeV (= 9.816 GeV) , (78) 
E bb (us; pQCDLamb.) = -0.355 ± 0.201 GeV , (m b = 4.382 GeV, k s / us = 0.306) ; (79) 

M r{is) (s; pQCDMS) = 9.060 + (1.183 - 0.193) + (0.399 - 0.263) + (0.330 - 0.439) + (0.475 - 0.475) GeV 

= 9.060 + 0.991 + 0.136 - 0.109 + 0.000 GeV (= 10.078 GeV) , (80) 
E bb (us; pQCDMS) = -0.617 ± 0.394 GeV , (m b = 4.505 GeV, k s/us = 0.248) ; (81) 



We can see from Eqs. (74)- (80) explicitly that for the chosen corresponding central values of the parameter k s / us the 
renormalon cancellation is exact in the last term (fifth term, named is(s)) in the sum for the soft mass Mx(is) (s), 
except in the case of MA/APT where k s / us = 1.0 was chosen and the cancellation in t 5 (s) is approximate. 



In the extracted results for m b , Eqs. (66 1- (73), the uncertainty from the ultrasof t (u s ) se ctor is rather strong, 
especially in the pQCD approaches. The us binding energy E bb (us; m b ; /z„ s ), cf. Eqs. (59)-(63), depends strongly on 



the s-us separation parameter k s / us , the central value of the latter being fixed by the renormalon cancellation of the 
soft mass Mx(is)(s) in the last term (or, in MA/APT, k s / us — 1). The renormalization scale /i us was not allowed 

to go below (/i us ) m in = 1.1 GeV, because at [i us s» 1.1 GeV the coupling As(ijl us ) (= A5(u us ) in our approach) 18 
in 2danQCD reaches a local maximum, indicating that the model may not necessarily give reliable values of A5 and 
A§ t \ below such scales. Therefore, in general, we will not use scales lower than 1.1 GeV, in 2danQCD and in pQCD, 



i.e., in the us binding energy expressions ( 60 )-( 61 ) . It turns out that the latter expression, Eq. (61), in 2danQCD 
reaches a local minimum at fi us slightly above 1.1 GeV. Therefore, we estimate the us part of the binding energy in 
2danQCD in the following way [we use the notation E bb (us; 

E bb (us; 2danQCD) = 1 [E bb (us; Q 2 S ) + (E bb (us)) min \ ± 1 [E bb (us- Q 2 S ) - (E bb (us)) min \ , (82) 



where the soft reference scale Q s was determined by the condition (55), and gave for the central values of the input 
parameters (sq and c-i of the second line of TablejlJ a s (M|;MS) = 0.1184; fn b = 4.353 GeV) the soft reference scale 
value Q s « 2.08 GeV. 



We recall that in 2danQCD we calculate in m q (mj l ; fj.^) the couplings An(fJ^) as sums of A. p (fi^) with p = n, n + 1, , 4; and in 

Eqq(s\m^; fj,^) and E q q(us\rn^; (i^s) the couplings «4 nj fc(/^) (and A. n ,h (/^D) as (derivatives of the) sums of A v (^) (and A-p(Hus)) with 
p = n, n + 1) • • • , 5; cf. discussion in Subsec. IIIBIand Sec. Ill II 
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In the two pQCD approaches, the us binding energy turns out to decrease monotonously when decreases below 
the soft reference scale Q s . P t of Eq. (54). For pQCD in the Lambert scheme (c2 = —4.76), with central values of the 
input parameters, the soft reference scale turns out to be Q s , P t ~ 1-73 GeV, and the us binding energy at fj, us = 1.1 
GeV reaches the value of —0.56 GeV; we use the following estimate for the us part of the binding energy: 



1 



1 



£ 65 ( US ;Lamb.) = - [E bb (us; Q 2 pt ) + E bb (us; l.l 2 GeV 2 )] ± - [E bb (us; Q 2 pt ) - E bb (us; l.l 2 GeV 2 )] 



(83) 



In pQCD in MS scheme, however, E bb (us; l.l 2 GeV 2 ) sa —1.5 GeV, which is exceedingly low and indicates failure of 
the method already at such scales, due to vicinity of Landau singularities of the running coupling. 19 Therefore, in 
MS we take as the minimal acceptable scale /j> us — 1.2 GeV, where E bb (us; 1.2 2 GeV 2 ) ~ —1.0 GeV (and Q s ^ pt ~ 1-87 
GeV) when the central values of the input parameters are used 



E bb (us;MS) 



1 



1 



'E bl {us- Ql pt ) + E bb {us- 1.2 2 GeV 2 )] ± - [E bb (us; Q 2 , pt ) - E bb (us; 1.2 2 GeV 2 )] 



(84) 



In the MA/APT case, on the other hand, Q s is fixed according to Eq. (57 1 and gives, for the central input parameter 
values, the value Q Sj (ma) ~ 1-60 GeV. In MA/APT no practical problems appear at scales fj, ua < 1.1 GeV. It turns 
out that E bb (us, MA/APT) reaches a moderate maximum value of 0.009 GeV at fi us ss 0.7 GeV for the chosen central 
value k s / us = 1.0. Therefore, we estimate the us part of the binding energy in MA/APT in the following way: 



E bb -(us; MA/APT) 



l r 



(E bb (us)) r 



E bb (us;Q 2 



s,(MA), 



± 



l r 



E bb (us;Q s ( MA) ) - [E bb (us)) 

1 



(85) 



In Eqs. (66l-(73), the uncertainties in m b originating from this determination of the us binding energy are denoted 
by the subscript (us). 

The described behavior of the us binding energy E bb (us; /i 2 s ) in the three models (2danQCD, and pQCD in the two 
schemes), as a function of the us renormalization scales fi us in the low-momentum regime, is presented in Fig. |3][a), 
and in the case of MA/APT in Fig. [^b). 




0.5 1.0 1.5 2.0 0.5 1.0 1.5 2.0 



FIG. 3: (a) The ultrasoft binding energy E bb -(us) (in GeV) as a function of the us renormalization scale fi U3 , for the central 
input values of 2danQCD model for the T(1S) system (mi, = 4.353 GeV; k s / us = 0.241; C2 = —4.76; so = 23.06): solid line 
for 2danQCD; dashed line for Lambert pQCD (c 2 = -4.76); dash-dotted line for the 4-loop MS pQCD (c 2 = 4.471). For all 
three curves the same nib and k s / us values are used. Included are also the properly rescaled 2danQCD couplings A$ and .4s, i 
as functions of the scale \i (= /i U3 ). (b) The same as in Figure (a), but for the MA/APT model (with c 2 = 4.471), with the 
corresponding central values of that model (mt, = 4.165 GeV and k s / ua = 1). In all the curves Nf = 3 was taken. 



On the other hand, there are related uncertainties for the extracted fn b , denoted by the subscript (s/us) in Eqs. ( 66 )- 
(|72| , which originate from the variation of the s-us separation parameter k s / us . The parameter k s / us was varied in 
such a way that the last (fifth, £s(s)) term in the series for the soft mass M T ( 1S )(s), cf. Eqs. ( 74 )-( 80 1, varies between 



It turns out that the Landau cut in the perturbative coupling a p t(/i 2 ) starts in the 4-loop (Nf = 3) MS scheme already at about 
fj, = Al ~ 0.607 GeV; while in the Lambert scheme [Eq. ( |10| with the central value c 2 = —4.76] it starts at a significantly lower scale 
fi = Al ~ 0.262 GeV. The 2danQCD coupling Ai(fi 2 ) has, of course, no Landau singularities, although it almost coincides with the 
Lambert pQCD coupling a p t(/J 2 ) at higher scales |^ 2 | > 1 GeV 2 . 
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the penultimate term t4(s) of these series, and its negative —ti(s), these two cases corresponding to the upper and 
the lower entry of (s/us) uncertainty of mj,, respectively. Only in the MA/APT case, where the exact renormalon 
cancellation was not achieved, the parameter k s / us was varied between 0.1 and 10, i.e., k s / us — l.Crtgg. 

The other uncertainty in the determination of m& comes from the uncertai nty of the A 2 scale. In MA/APT it comes 
from A 5 = 0.260 ± 0.030 GeV and is denoted by the subscript (A 2 ) in Eq. g. In 2danQCD and in the two pQCD 
approaches (Lambert scheme and MS scheme), this uncertainty comes from a s (M^;MS) = 0.1184 ±0.007 [44< and is 
denoted by the subscript (a s ) in Eqs. (68), (70), (72). 
Yet another uncertainty of rhb 



in 2danQCD and in Lambert schemepQCD, comes from the va riati on of the 
(Lambert) renormalization scheme parameter C2 = — 4.76^0 97! c f- Table I 20 and is denoted in Eqs. (68) and (70) 
by the subscript (02). On the other hand, the scheme in MA/APT was fixed by the underlying pQLT) solution, 
Eqs. |l0|-|l3j), with c 2 = c 2 (MS, Nf = 3) « +4.471 (and c 3 = 4 _1 /^~ 2 for 3 ^ 3 )- In (4-loop) MS pQCD, the 

-4.471; c 3 = c 3 (MS, N f = 3) « +20.99; c 3 = 



scheme parameters are fixed in the usual way: c 2 = c 2 (MS, Nf = 3) 
for j > 4. 

When we vary the soft renormalization scale /z 2 around the central value (Q 2 ) CO ntr. of the soft reference scale 
[determined by Eqs. (54), (55) and (57) in pQCD, 2danQCD and MA/APT, respectively, for central values of the 
input parameters M b , k s / us , etc.], so that ^ 2 = 2(<3 2 ) contr . or fi 2 s = (l/2)(Q 2 ) ccntr 
variations of extracted m& as denoted by the subscript (^i 2 ) in Eqs. (66)-(72). 



this results in the corresponding 



Finally, the uncertainty Smb(m c + 1 0) = ±0.005 GeV due to nonzero c mass, Eq. (140 1 
in the uncertainties +0.005 GeV of TUb, denoted by the subscript (m c ) in Eqs. ( 66 )-( 72 ) 
We see in Eqs. (pl-d^f 



cf. also Eqs. (46 )-(49 )] results 



that the largest resulting uncertainty in the determination of rUb is the one originating 
from the uncertainty of the determination of the us binding energy (except in MA/APT where l-E^its)] values 
are small). These uncertainties are larger in the two pQCD approaches, due to the influence of the (unphysical) 
Landau singularities in the running couplings. Furthermore, the effect of the us sector in 2danQCD and in the pQCD 
approaches was to increase the predicted value of mi, because the us binding energies are in these cases significant and 
ncg ative, cf. Fig.^a). If we had ignored the existence of the us sector, i.e., if we had used in the entire binding energy 
E b b simply a common soft renormalization scale fi ~ Q s , the predicted values of rf%b in 2danQCD and in Lambert and 
MS pQCD approaches would have decreased, by —0.068, —0.091, and —0.177 GeV, respectively, as can be deduced 
from the its-origin uncertainties in Eqs 



(68) 



72). On the other hand, in MA/ APT the choice fj, us = Q 



s,(MA) 



would 



}y 0.002 GeV, cf. Eq. (66), basically because the values of \E bb -(us)\ in 



only slightly increase the central value of mj 
MA/APT are much smaller, cf. Figs.^a),(b) 

For better visibility, we present the results for the central extracted values of WTb of the aforementioned four methods 



in Table nil and for various uncertainties drrib in Table III 



TABLE II: The extracted central values of mb in the four models, for the central input parameter values (with the total 
uncertainties Sffib). Included are the corresponding input parameter k s / us , and the resulting scales: soft reference scale Q s 
(soft ren orm a lizat ion scale is taken /j, s = Q s ); soft mass Mxns)(s); averaged ultrasoft energy E bb -(us) and its ambiguity 8E bb -(us) 
[cf. Eqs. ( 82 1- ( 85 1] . All the scales are given in GeV. Note that M T ( 1S )(s) + E bb -(us) = 9.460 GeV, i.e., the physical mass M r (is)- 



model 


mb(5m,b) 


h , 

"'S /US 


Q s 




E bb(us) 


SE bi 


MA/ APT 


4.165(±0.022) 


1.000 


1.599 


9.454 


0.006 


T0.003 


2danQCD 


4.353(±0.084) 


0.241 


2.084 


9.814 


-0.353 


±0.151 


pQCD Lamb. 


4.382(±0.111) 


0.306 


1.729 


9.816 


-0.355 


±0.201 


pQCD MS 


4.505(±0.231) 


0.248 


1.869 


10.078 


-0.617 


±0.394 



We wish to address here briefly the question of nonperturbative (NP, higher-twist) contribution to the quarkonium 
mass. In the heavy quark system such as bb, the nonperturbative (NP) contribution can be estimated, and in the 
leading order it comes from the gluon condensate and is given by |78j 



E bl {us) 



< NP > « m 67 r 2 



624 /4tt 
425 



1 



a pt(ML 



(86) 



We recall that in 2danQCD all the other dimensionless parameters vary with C2 according to Table |7J the four relations Sj(so) an d 
fj(so) (for j = 1, 2) are fixed from the (four) conditions Ai(Q 2 ; C2) — a pt (Q 2 ; C2) ~ (A 2 /Q 2 ) 5 for \Q 2 \ > A 2 ; the scale A is fixed by the 
condition q s ( M|; M S) = 0.1184; and the parameter so (= M^/A 2 ) is then fixed by the condition r T (V + A) = 0.203. For details, we 



refer to Subsec. II B Higher Lambert scheme parameters are fixed by: Cj 



"Vcj" 2 for j > 3, cf. Eqs. |lu|-|l3| 
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TABLE III: The uncertainties Sfrib of the extracted value of frib coming from various sources: (1) from the evaluation of the us 
sector; (2) from the variation of the s-us separation parameter k s / us ; (3) from the variation of a s (or, in MA/APT: variation 
of A); (4) from the variation of the C2 parameter (in 2danQCD, and in pQCD in the Lambert scheme); (4) from the variation 
of the soft renormalizatio scale u s ; (5) from the uncertainty of (<5m(,)m c of Eq. (40 1. See the text for details. 



model 


5irib(us) 


5m,b(s i 'us) 


{J^s/us ) 


Sm b {a s ) 


5m b (c 2 ) 


(<*) 


Sm b (u 2 s ) 


5m b {m c ) 


MA/ 


+0.002 


+0.005 


(1.0+9.0) 


-0.019 




(") 


-0.005 


-0.005 


APT 


-0.002 


-0.004 


(1.0-0.9) 


+0.020 




(") 


+0.003 


+0.005 


2dan- 


-0.068 


+0.015 


(0.241-0.100) 


-0.005 


-0.023 


(-4.76+2.66) 


+0.017 


-0.005 


QCD 


+0.071 


-0.016 


(0.241+0.173) 


+0.005 


+0.034 


(-4.76-0.97) 


-0.025 


+0.005 


pQCD 


-0.091 


-0.013 


(0.306+0.194) 


+0.010 


+0.027 


(-4.76+2.66) 


-0.002 


-0.005 


Lamb. 


+0.097 


+0.017 


(0.306-0.119) 


-0.010 


-0.008 


(-4.76-0.97) 


-0.041 


+0.005 


pQCD 


-0.177 


-0.082 


(0.248+0.642) 


+0.031 




(") 


-0.004 


-0.005 


MS 


+0.200 


+0.084 


(0.248-0.181) 


-0.027 




(") 


-0.075 


+0.005 



The factor l/a£ t (^ a ) in (4-loop) MS pQCD is unreliable for realistic us scales n us s» 1 GeV, due to the vicinity of the 



Landau singularities (cf. footnote 19). In 2danQCD, for purposes of estimation, we replace l/a 4 t (^ 2 s ) by l/A^i^s) 

or by l/Aids^s)- In the interval (1.1 GeV < u, us < 1.3 GeV) we have A4 ~ A4 ~ 10~ 4 ; the couplings A4 and A4 
cover in this interval the values between 0.5 x 10 -4 and 2 x 10~ 4 . For these values, and using the central value of the 
gluon condensate ((a s /7r)G 2 ) = 0.009 GeV 2 [75] (cf. also Refs. [SO] and [53]), and for M b = 4.3 GeV, we obtain the 
following estimate: 



E bl (us)^=0M±™l GeV 



(87) 



This effect is relatively small and has large uncertainties. If we take it into account, then the central extracted values 
of Mb in this Subsection decrease somewhat, the decrease being (Am^Wp) s» (— 0.025^o.'qio) GeV. 

We wish to dwell also on the following aspect: as seen from the results of this Subsection, the extracted values 
nib and various uncertainties Sfrib are similar in 2danQCD and the corresponding Lambert pQCD. The main reason 
for this is the fact that the scheme parameters (c,-, j > 2) are the same in both frameworks, and that the two 
corresponding running couplings practically merge at high momenta \Q 2 \ > A 2 : Ai(Q 2 ) — a pt (Q 2 ) ~ (A 2 /Q 2 ) 5 . On 
the other hand, the method of evaluation of the binding energies E qq (s) and E qq (us) and of the pole mass m q was 
different in the two approaches. In (Lambert )_pQCD, these quantities were evaluated as truncated power series in 
powers of a p t [the latter determined by Eq. (10)], and the terms et pt In a pt in E qq were also evaluated directly by 
using the corresponding value of a p t - cf. Eqs~147) and (46). On the other hand, the model 2danQCD is a case of 
general analytic QCD, with the basic coupling A~i(Q 2 ) given in Eq. (18), and the evaluations of the higher power 
analogs (a pt )an = A n and of the logarithmic term analogs (a pt ln fc 



a pt J 



A n ,k are more complicated in such a case, 



because linear combinations of the logarithmic derivatives A v are needed for such a construction [A v for integer v 
is given in Eq. (27), and for noninteger v in Eq. (31)] - cf. Eqs. ( 28 )-( 30 1 and (34). These latter relations involve 
truncations. For terms in E qq , Eqs. (50) and (61|-(64), we truncated the mentioned relations ( 28 )-( 30 1 and (34) 
at A5 (~ dp t ) because perturbation expansion of E qq is known up to ~ et pt terms; for terms in the pole mass m q , 
Eq. (53), we truncated the relations (28)-(30) at A4 (~ a pt ) because perturbation expansion of m q is known up to 
~ a pt terms [when the estimated r$ term, Eq. (45), is included]. We expect that the difference between the results m& 
(and Sfrib) in the 2danQCD and in Lambert scheme pQCD would decrease even more if we used in both frameworks 
the same method of evaluation. For example, in (Lambert) pQCD we could use for the powers of a pt the linear 
combinations Eqs. ( 24 )-( 26 1 in terms of the logarithmic derivatives a p t,n, Eq. (20), where the linear combinations 
Eqs. (24)-(26) would be truncated at a p t,5 for the powers in E qq , and at o p t,4 for powers in m q . We note that the 



logaririmiicderivatives in Lambert pQCD and 2danQCD differ at large momenta (|Q 2 | > A 2 ) also negligibly; namely, 
if k is integer, we can apply k — 1 times the derivatives Q 2 dl(dQ 2 ) to the relation A\(Q 2 ) — a P t(Q 2 ) ~ (A 2 /Q 2 ) 5 
and obtain A k {Q 2 ) - a pt . k (Q 2 ) - (A 2 /Q 2 ) 5 (cf. also Figs.hM There is a nontrivial technical problem, though, in 
the evaluation of the logarithmic terms a™ t In a pt in E qq along these lines in (Lambert) pQCD - for this, we would 
need, according to Eq. (34), to evaluate derivatives of the type d k ~ s a p t. u + m /di ,k ~ s at v = n integer, but this implies 
that 



is needed also for v noninteger (since variation of v around n is needed), implying that we would need 



to know the logarithmic derivatives H v t,v at noninteger v, i.e., a relation of the type (31 ) would be needed for pQCD. 
This can be probably achieved by using in the integral ( |3l"| ) the pQCD discontinuity function for pi(cr) and extending 
the integration down to a — — A^ — 0, i.e., down to the beginning of the (unphysical) Landau cut, in the spirit of 
Eq. |l]). In the present work, we did not perform this calculation for (Lambert scheme) pQCD, though. 
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B. Charm mass extraction 



In this case, qq = cc, the quarkonium mass is now Mj/^^g) = 3.0969 GeV [44]. We basically repeat the analysis 
as in the case of bb. There are some differences, though: 



• The relation (41) is now at Nf — 3; therefore, the transition to the relation at the soft renormalization scales, 
Eq. (46), which is also at Nf = 3, has now no threshold transition complication. 



The typical soft reference scales [Eqs. (54) and (56)] are now significantly lower, Q s w 1 GeV or even lower 
(in the bb case we had: Q s « 2 GeV). This, in conjunction with our suggestion that the considered models 
2danQCD and pQCD (in the Lambert and MS schemes) are not necessarily to be trusted at scales below 1.1 
GeV, implies that the typical soft renormalization scales /z s should be chosen significantly higher than Q s . We 
choose /ig = (3 ± 1)Q 2 in these three models. In this case, the lowest possible soft scale, (^ s ) m i n = V^Q S 
is slightly above 1.1 GeV. In MA/ APT, the soft renormalization scale will also be varied in this way, giving, 
however, somewhat lower central value for the renormalization scale, fi s — V^Qs^ma) ~ 1-004 GeV. 21 In general, 
by keeping now the soft renormalization scale /i s higher than the (low) soft reference scale Q s , the uncertainties 
(<5m c )( M 2) of the extracted mass m c due to the scale variation /x 2 = (3 ± 1)<3 2 may be underestimated. 

In general, the cancellation of the leading renormalon now implies for the s-us separation parameter k s / us such 
values for which the absolute values of the us binding energy E c5 (us) are significantly smaller than in the bb 
case, and consequently, the us ambiguities are smaller. Since we now use for the central choice of the soft 
renormalization scale u 2 s — 3Q 2 , the us binding energy is calculated in 2danQCD and in pQCD in the following 
way: 

E c5 (us) = 1 [E c5 (us; 3Q 2 ) + E c5 (us; l.l 2 GeV 2 )] ± l - [E c5 (us; 3Q 2 ) - E c5 (us; l.l 2 GeV 2 )] , 



where in 2danQCD the soft reference scale Q s is determined via Eq. (56) and in pQCD Q s = Q s , P t 



determined via Eq. (54). In MA/APT, we do not appear to have practical problems at low scales fj, < 1.1 GeV, 
and the us energy as function of low scales u us turns out to have a moderate local maximum and a moderate 
local minimum; hence we use the formula 

E cS (us; MA/APT) = i [E C c(us) max + E c£ (us) min ] ± i [£^(us) max - E cS (us) min ] , (89) 

where these values are quite small: in the central case (fn c = 1.262, k s / us = 1.0), the local maximum is 
E C c{us) max « +0.003 GeV and is reached at u us ss 0.69 GeV, and the local minimum is E C c(us) m i n ~ —0.005 
GeV and is reached at very low scale u us w 0.10 GeV. 

• The exact renormalon cancellation requirement in MA/ APT now does not give an exceedingly small value of 
the s-us separation parameter k s / us , but it it still quite small, k s / us ~ 0.02. We decide, for presentation reasons, 
to vary in MA/APT the parameter k s / us again around its central chosen value 1.0, in the interval between 0.1 
and 10., just as it was done in the bb case of MA/APT. 

The described behavior of the E c5 (us; yLt 2 s ) in the analytic 2danQCD model and in pQCD in the two schemes, as a 
function of the ulttrasoft scales fi us , is presented in Fig. |4ja), and in the case of MA/APT in Fig. |4jb). Comparing 
Figs. ^a) and (b), we see that \E c5 (us) \ - lO^-lO" 1 GeV in 2danQCD, and - 10~ 3 GeV in MA/APT. Furthermore, 
comparing with the corresponding curves in Figs. [3^ a) and (b) for the bb, we can see that in 2danQCD and in the two 
pQCD approaches, the absolute values |-E ce (us)| are by almost two orders of magnitude smaller than \E b j } (us)\, prin- 
cipally because the renormalon cancellation gives us in the two cases significantly different s-us separation parameter 
values: k s / us {cc) = 4.00 and k s i us (bb) — 0.241. Further, in general the values of \E q q(us)\ in pQCD are larger than 
in 2danQCD, this being true even more so in MS pQCD than in Lambert scheme pQCD, basically because of the 
relatively high branch cut scale A 2 at which the Landau singularities cut start in the MS case (cf. footnote 19 1 . In 
the MA/APT case, the values of \E q q{us)\ are quite small in both cases (q — c, b), being in cc case smaller by almost 
one order of magnitude. All this is reflected also in the numerical results of this and of the previous Subsection. 



If we used in MA/APT a lower definition of the central renormalization scale, fi s = Q s ,(mA) (~ 0.58 GeV), the predicted central value 
of m c would go up by only 0.005 GeV. 

Note that in 2danQCD in the cc case the condition (1556 cannot be fulfilled. 
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FIG. 4: (a) The ultrasoft binding energy E C c(us) (in GeV) as a function of the us renormalization scale fi us , for the central 
input values of 2danQCD model for the J/ip(lS) system (m c = 1,266 GeV; k s / us — 4.00; C2 = —4.76; so = 23.06): solid 
line for 2danQCD; dashed line for Lambert pQCD (c 2 = -4.76); dash-dotted line for the 4-loop MS pQCD (c 2 = 4.471). For 
all three curves the same m c and k a / us values are used. Included are the properly rescaled 2danQCD couplings AsifJ-ts) an d 
A5 : i((J.us) as functions of the scale fi us . (b) The same as in Figure (a), but for the MA/APT model (with c 2 = 4.471), with 
the corresponding central values of that model (m c = 1.262 GeV and k s / us = 1). In all the curves Nf = 3 was taken. 



The resulting extracted values of m c are 

m c (MA/APT) = {1.262 + 0.002 (us) + 0.002 (s/ „ s) + 0.011 (A 2) + 0.002^2)} GeV 



1.262 ± 0.012 GeV , with : \/3Q S:(M a) = 1-00 GeV , k s/us = 1.0 



(90) 
(91) 



m c (2danQCD) = i 1.266 ± 0.003 (us) ± 0.007 (s /„ s) T 0.005 (c 



-0.014 
-0.003 



(c 2 ) 



-0.002 
-0.005 



(Ms) 



1.266 ± 0.017 GeV , with : V3Q S = 1.42 GeV , k s/us = 4.0 . 



m c (pQCDLamb. 



1.265 ±0.001 



(us) 



+0.021 
-0.021 



(s / us) 



T0.004 (Qa) ±0.000 (c2) 



-0.003 
-0.015 



1.265 ± 0.027 GeV , with : VSQ S = 1-38 GeV , k s/us = 5.59 . 



GeV (92) 
(93) 

GeV (94) 
(95) 



to c ( P QCDMS) = i 1.272 t0.011 ( „ s) + J t 0.002 (Qs) + 1 ^ GeV (96) 

= 1.272 ± 0.078 GeV , with : V3Q S = 1.58 GeV , k s/us = 3.08 . (97) 

In order to see more clearly the renormalon cancellation, we present below, as in Subsec. |IV A[ the sum for the soft 
mass A^j/^(is)(s), combining in each parenthesis the (positive) term ~ a" from 2m c (m^; /ijf) and the (negative) term 
~ a n+1 from the soft binding energy E cS (s; m 2 c ] /i^) (n — 0,1, ... , 4), for the central input values of parameters m c , 
Ms ( = V%Qs) and k s / us ; separately we present below also E c5 (us). 

2.524 + (0.412 - 0.056) + (0.227 - 0.081) + (0.125 - 0.071) + (0.045 - 0.026) GeV 
2.524 + 0.356 + 0.146 + 0.054 + 0.020 GeV (= 3.099 GeV) , (98) 
-0.001 ± 0.004 GeV , (m c = 1.262 GeV, k s/us = 1.0) ; (99) 



M J/q p ilS )(s; MA/APT) 

E cS (us; MA/APT) 
M J/v , (ls) (s;2danQCD) 

E cS (us; 2danQCD) 
M.j/iP(iS){s; pQCDLamb. 
E cS (us; pQCDLamb. 



2.532 + (0.349 - 0.053) + (0.328 - 0.171) + (0.253 - 0.150) + (0.979 - 0.979) GeV 
2.532 + 0.296 + 0.157 + 0.103 + 0.000 GeV (= 3.088 GeV) , (100) 
0.009 T 0.006 GeV , (m c = 1.266 GeV, k s/us = 4.0) ; (101) 

= 2.530 + (0.354 - 0.061) + (0.303 - 0.141) + (0.404 - 0.305) + (0.662 - 0.662) GeV 
= 2.530 + 0.293 + 0.161 + 0.099 + 0.000 GeV (= 3.083 GeV) , (102) 
= 0.013 T 0.003 GeV , (m c = 1.265 GeV, k s/us = 5.59) ; (103) 
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TABLE IV: The extracted central values of m c in the four models, for the central input parameter values (with the total 
uncertainties Sm c ). Included are the corresponding input parameter k s / ua , and the resulting scales: soft renormalization scale 
Ms — V3Q S ; soft mass Mj/^^s){ s ): averaged ultasoft energy E C c(us) and its ambiguity 8E C c(us) [cf. Eqs. f 88 )-( 89 1] . All the 
scales are given in GeV. Note that Mj/^(!S)(s) + E C c(us) — 3.097 GeV, i.e., the physical mass Mj/^^s)- 



model 


m c (Sm c ) 


h , 

<"3 1 US 






E C c(us) 


5£ C c 


MA/ APT 


1.262(±0.012) 


1.00 


1.004 


3.099 


-0.001 


±0.004 


2danQCD 


1.266(±0.017) 


4.00 


1.423 


3.088 


+0.009 


+0.006 


pQCD Lamb. 


1.265(±0.027) 


5.59 


1.382 


3.083 


+0.013 


+0.003 


pQCD MS 


1.272(±0.078) 


3.08 


1.585 


3.131 


-0.034 


±0.024 



TABLE V: The uncertainties 5m c of the extracted value of m c coming from various sources: (1) from the evaluation of the us 
sector; (2) from the variation of the s-us separation parameter fc s / us ; (3) from the variation of a s (or, in MA/APT: variation 
of A); (4) from the variation of the C2 parameter (in 2danQCD, and in pQCD in the Lambert scheme); (4) from the variation 
of the soft renormalizatio scale /i s ; See the text for details. 



model 


Sm c {us) 


Sm c (s/us) 


(^S / liS ) 


Sm c (a s ) 


<5m c (c 2 ) 


(ca) 


<$m c (/4) 


MA/ 


-0.002 


-0.002 


(1.0+9.0) 


-0.011 




H 


-0.002 


APT 


+0.002 


+0.002 


(1.0-0.9) 


+0.011 




(-) 


+0.002 


2dan- 


+0.003 


+0.007 


(4.00-3.61) 


-0.005 


-0.014 


(-4.76+2.66) 


-0.002 


QCD 


-0.003 


-0.007 


(4.00+39.5) 


+0.005 


+0.003 


(-4.76-0.97) 


+0.005 


pQCD 


+0.001 


+0.021 


(5.59-5.355) 


-0.004 


0.000 


(-4.76+2.66) 


-0.003 


Lamb. 


-0.001 


-0.022 


(5.59+148.4) 


+0.004 


0.000 


(-4.76-0.97) 


+0.015 


pQCD 


-0.011 


+0.066 


(3.08-2.86) 


-0.002 




H 


-0.003 


MS 


+0.011 


-0.075 


(3.08+66.9) 


+0.002 




(") 


+0.017 



Mj /v , ( i S) (s;pQCDMS) = 2.544 + (0.368 - 0.066) + (0.348 - 0.163) + (0.455 - 0.356) + (0.778 - 0.778) GeV 

= 2.544 + 0.302 + 0.185 + 0.099 + 0.000 GeV (= 3.131 GeV) , (104) 
E c£ (us; pQCDMS) = -0.034 ±0.024 GeV , (m c = 1.272 GeV, k s / us = 3.08) ; (105) 

As in the previous Subsection in the case of bb, we present now for the case cc, for better visibility, the results for 
the central extracted values of fn c of the four models in Table |IV| and for various uncer tain ties 5m c in Table |V} 



The nonperturbative (NP) contribution coming from the gluon condensate, cf. Eq. (86) for the bb system in the 
previous Subsection, is for the lighter cc system unreliable, since the next-to-leading corrections are in this case large 
and tend to make the result unreliable, cf. Ref. [8"T1 . 



V. SUMMARY 



We evaluated, in two analytic QCD models and in perturbative QCD (pQCD, in two schemes), the quark-antiquark 
binding energies (up to N 3 LO) and masses of ground state qq quarkonia (q = b,c), as functions of the quark mass m q 
(= m q (u 2 = to 2 )), also called MS quark mass. In analytic QCD models the QCD running coupling A\(Q 2 ) has no 
unphysical (Landau) singularities in the Q 2 -plane. One analytic model used was the Minimal Analytic model (referred 
to also as: Analytic Perturbation Theory, MA/APT) of Shirkov, Solovtsov, Solovtsova, and Milton, Refs. [TTHHj. 
The (analytic) QCD running coupling in this model, A < l AA \Q 2 ), differs nonnegligibly from the corresponding (same 
scheme) pQCD coupling a pt (Q 2 ) = a s (Q 2 )/n even at large momenta \Q 2 \ > A 2 : y^ MA) (Q 2 ) - a pt (Q 2 ) ~ (A7Q 2 ) 1 . 
The other analytic QCD model used was the so called 2-delta analytic QCD (2danQCD), Ref. [35]; which basically 
coincides with pQCD in the same (Lambert) scheme at high momenta \Q 2 \ > A 2 : A\(Q 2 ) — a p t(Q 2 ) ~ (A 2 /Q 2 ) 5 . 
The use of the analytic QCD models was motivated by the fact that the typical soft (s) momentum scales Q s in the 
ground bound states of quarkonia are low (Q s ps 2 GeV and 1 GeV, for bb and cc, respectively), and that the typical 
ultrasoft (us) momentum scales Q us are even lower. This, in conjunction with the fact that Landau singularities of 
the pQCD coupling a p t(Q 2 ) reach relatively high momenta: Q f=s 0.61 GeV in the usual (4-loop) MS scheme (with 
c 2 = /3 2 //3 = 4.471), and Q « 0.26 GeV in the Lambert scheme (c 2 = -4.76). 

In the course of the calculations, we included a renormalon-based estimate of the N 3 LO term (~ a pt ) of the 
quark pole mass m q in terms of a pt and m q , cf. Appendix This allowed us to fix the s-us separation of the 
binding energy, E qq = E qq (s) + E qq (us), via the requirement of cancellation of the leading infrared renormalon in 
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the sum M qq (s) = 2m q + E qq (s). In the latter sum we used as the input parameter the mass m q which is free of 
the renormalon ambiguity, and the coefficients of the perturbation expansion in a p t(Ms) a ^ the same common soft 



pt 



renormalization scale fj,j., where for the quark pole mass and for the soft binding energy terms up to ~ a* t and 
were used, respectively. It turned out that in MA/APT the renormalon cancellation condition cannot be exactly, 
but only approximately, fulfilled. In the case of the bb system (T(ls)), this condition gave in 2danQCD and in 



pQCD the values of the s-us separation parameter k s / us , Eq. (59), between 0.24 and 0.31. This led to significant 



(and negative) values of the us part E b i{us) of the binding energy, cf. Fig. [3](a). In pQCD these values had a 
particularly large uncertainties due to the vicinity of Landau singularities. These effects of the us regime led, in 
comparison with the usual perturbative approaches in the literature, to higher central values of the bottom quark 
mass: m& w 4.35,4.38,4.51 GeV in 2danQCD, Lambert p QC D, MS pQCD, respectively. Further, these effects led to 
large uncertainties, especially in MS pQCD - cf. Tables |ll| and III If we used in E b i{us) the same (soft) renormalization 
scale as used in E bb (s) (u ~ Q s ~ a s ni b ), the extracted central values of m& would decrease, by about —0.07, —0.09, 
—0.18 GeV, respectively. The 2danQCD analytic model gave us, among other things, clues about how to estimate 
the large negative values of E bb (us) and the uncertainties thereof, not just in 2danQCD, but also in pQCD in both 
schemes. The other, less significant, uncertainties in the extraction of the value of m b , come from the uncertainty of 
the s-us separation parameter k s / us , from the uncertainty of the world average of the QCD coupling a s (M|;MS), 
from the renormalization scheme and renormalization (soft) scale variations (ci and u s ), and from the uncertainty of 
the finite charm quark mass effect at fj, s scales. On the other hand, MA/APT has in general quite suppressed couplings 
in the infrared, so that the us sector does not affect it significantly. Our results for m b are given in Eqs. (66)- (73 1, 
and in Tables [ED and HU 



In the cc system (J/ip(lS)), the analysis was similar. Surprisingly, the infrared renormalon cancellation condition 
gave in this case for the s-us separation parameter k s / us significantly higher values, between 3.1 and 5.6 in 2danQCD 
and in pQCD (in MA/ APT the cancellation condition was again not achieved). As a result, despite having in general 
even lower typical us scales fi us than in the bb case, the us parts E cS (us) of the binding energies have in general 
significantly smaller absolute values than in the bb case, cf. Figs, ga) and^a). This resulted in smaller influence of 
the us sector on the extraction of m c , and also in smaller uncertainties 5m c from this sector, than in the case of bb 
(Mb and Sfrib) - cf. Eqs. (90)-(97), and Tables IV and[v} The obtained central extracted values of m c w 1.27 GeV are 
thus close to those in the literature. In MA/ APT the results for m c are also similar. 

The main conclusion of this work appears to be the following: analytic QCD approaches which at high energies 
follow the pQCD behavior closely (such as 2danQCD analytic QCD) indicate that the ultrasoft regime in the bb 
system is important. Our approach indicates that the influence of the ultrasoft regime is underestimated in the usual 
pQCD approaches (especially in the MS scheme) which consider the soft and ultrasoft energies to be roughly equal 
and consequently give lower values of the extracted masses nib (nib < 4.3 GeV). 
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Appendix A: Ground state quark-ant iquark binding energy 

According to Refs. [71 1581 164] . the perturbation expansion of the quark-antiquark binding energy E qq , in terms of 
the quark pole mass m q and a s (u 2 ) = a s , is 



E qq = E? + SE[ 1] + 8Ef> + 5E[ 
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+21n[2]) + — C A C F N f T F + C 2 F N f T F ( 
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The following notations were used: 
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The RGE coefficients /?,■ are in MS scheme 
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The constants ai and are 
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The value for the constant 0,3 asociated with the three-loop soft contribution was obtained in Refs. [3] and is 
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18 - 67V, 2 , + 7V£ 
96iV| 



(2) 

a < 3 > = - 



C F T F 



14002 416C(3)\ 



81 



C A T\ 



( 12541 64tt 4 368C(3)\ 



V 243 



135 



8000T| 
729 



(A13) 



Appendix B: Renormalon-based estimate of ~ ait coefficient 



The term r% in the expansion of m q /rn q in Eq. (41 1 can be estimated by a method closely related with the approach 
presented in Sec. II of Ref. [74]. The pQCD version of the sum in Eq. (41 ) can be reexpressed in terms of a pt (ii 2 ) at 
any other renormalization scale /i 2 



Tflq o 

ri(fi 2 ) = «i + l3 L m (n 2 ) , 

r 2 (/i 2 ) = k 2 + (2«i/?o + /?i)i m (^ 2 ) + /3o^m(M 2 ) , 

(4/3)«x = 6.248/3 - 3.739 , 

(4/3)k 2 = 23.497/3^ + 6.248/?i + 1.019/3 - 29.94 , 



(Bla) 

(Bib) 
(Blc) 
(Bid) 
(Ble) 



where L rn (fi ) = ln(/i /m q ), while fio(Nf) and /3\(Nf) are the renormalization scheme independent coefficients (given 



just after Eq. ( 13 >). Here, 2V/ = Ng is the number of light active flavors (quarks with masses lighter than 



25 



Since r\ and r 2 are explicitly known, the Borel transform Bs(b) is known to order ~6 2 



ri{n 2 ) r 2 (/x 2 ) 

nj 0+ 2\pi 



0(b 3 ) 



(B2) 



The function Bs(b) has renormalon singularities at b = 1/2,3/2,2,...,— 1,-2,... [82-84]. The behavior of Bs near 
the leading infrared (IR) renormalon singularity b — 1/2 is determined by the resulting renormalon ambiguity of m q . 
This ambiguity 5m q is a (QCD) scale which, having the dimension of energy and being renormalization scale and 
scheme indepenent, must be proportional to the QCD scale Aqcd^ 5m q — const x Aqcd [85] . This scale can be 
expressed in terms of a p t(n 2 ) and \i (fi being any renormalization scale) in the form 



where 



Cl 



A = const x fiexp 



2/3 ~ 2/3 2 ' 
(c 2 ~c 2 ) 
{2f3 fv ' 



1 

2/?oa P t(M) 



1 + J2 (20o) k "{v ~ 1) • • • (" - * + l)cfcflpt(M) 



fe=i 



(B3) 



C3 



6(2/3o) 6 ^(^-l)(^-2 



(B4a) 

2(2/3b)M,-l) ^ - C2)2 - ^ - 2ClC2 + C3) ] ' (B4b) 
r M - c 2 ) 3 - 6/3 (c 2 - c 2 ){c\ - 2 Cl c 2 + c 3 ) + 8/? 2 (c 4 - 3c 2 c 2 + c\ + 2 Cl c 3 - c 4 )](B4c) 



The above constants, expressed in terms of /3o and of Cj = j3j/ (3q, appear in the expansion of the residue of the Borel 
transform Bs(b) at the pole b = 1/2 



B s {b;n) = N m ir 



m q (1 - 2b) 1+v 



l + £c fe (l-26)* 



k=l 



B 



(an.) 



(&;m) > 



(B5) 



where -Bg" 1 '^ (6; /x) is analytic on the disk |6| < 1 and can be expanded in powers of b We can assume that the 
coefficients Cfe are known up to k — 3, because the coefficient C4 = Pi/Po ( m MS scheme) is known to a large degree 
by Pade-related methods of Ref. (5(5] 



1 



Pa = ^(A 4 + BiN f + dNj + DiNf + E 4 Nj) 



(B6) 



with A± = 7.59 x 10 5 , B A = -2.19 x 10 5 , C 4 = 2.05 x 10 4 , L> 4 = -49.8, and £ 4 = -1.84. This gives c 4 = 123.7 for 
Nf — 3, c 4 = 97.2 for iV/ = 4, and c 4 = 86.2 for iV/ = 5. The residue parameter N m can be determined with high 
precision by using the idea of Refs. [ST], i.e., by calculating (cf. Refs. [7211751155] ): 



where 



1 it 1 

N m = r ^-Rs(b = 1/2) , 

[I 7T 



i? 5 (&;M) = (i-26) 1+ ^ s (&;^) , 



(B7) 



(B8) 



and the first coefficients in the expansion in powers of b of this quantity are known from the known coefficients T\ and 
r 2 . We can use a combination of truncated perturbation series and Pade approximants [1/1] for R$(b), as presented 
in Ref. [55] , and obtain 



0.575(iV/ = 3) 



0.555(iV) = 4) , 0.533(iY/ = 5) 



(B9) 



with the uncertainties in th ese v alues of roughly ±0.020. 

In the bilocal expansion (B5), the analytic part Bg (b;p,) can be taken as a polynomial in b, i.e., a truncated 
expansion in powers of b. The coefficients of the latter expansion can be related with rj(u 2 )'s by equating the 
expansion of Eq. ( B5 ) in powers of b with the expansion ( B2 ) , resulting in 



B 



(an.) 



(6; a*) 



fe>i h - p o 



3 m « ^0 



T(u + k + l- n) 

'■' ro + i-n) 



(BlOa) 
(BlOb) 



2G 



where, by convention, r$ — cq — 1. The numbers c n of Eqs. ( |B4| , which enter the sum in Eq. ( B10b[ ), are known only 
up to n = 3, because, in MS, only up to k — 4 are reasonably known (04 approximately, as mentioned). For Nf — 3, 
these values are: c x = -0.1638, c 2 = 0.2372, c 3 = -0.1205 (and: v = 0.3951). For N f = 4, they are: c x = -0.1054, 
c 2 = 0.2736, c 3 = -0.1610 (and: v = 0.3696 ). And for N f = 5 they are: c x = 0.0238, c 2 = 0.3265, H 3 = -0.2681 (and: 
v = 0.3289). Therefore, the sums in (BlOb) are truncated at n = 3. 



Theoretically, the pole closest to the origin in Bg '(b\n) is at b = —1, at least in the large-/3o approximation 



2:i 



Since in B^'Q); ju,) the coefficients h^\n 2 ) for k = 0, 1,2 are known (because r x (fi 2 ) and r 2 {^ 2 ) are known), we can 
construct the Pade [1/1] ( a „o (&) and check the pole of it. It turns out that this Pade, at the natural scale fx = m„, has 

-0.96, 



the pole at b 



-1.09, 



-1.12, for N 



f 



[2/1] 



i (6) have the pole at b 



3,4,5, respectively, reflecting correctly the theoretical expectation. 24 

m q ) by requiring that the Pade 
1. This gives us 



We can extend further this reasoning and obtain the next coefficient (at /i 



4 m) (m?) = -25.18(7V / 



3) , -28.28(iV/ = 4) , -35.62(7V / = 5) 



If constructing with these values of the other possible Pade approximant of index 3, namely [1/2] , 



(Bll) 
i (6) , it turns 

out that the nearest to origin pole of such Pade is th en at b = —1.003, —1.001, —1.008, for Nf — 3, 4, 5, respectively. 

This indi cates t hat the obtained values of Eq. (Bll), are consistent. 25 Using these values, we obtain from the 

relations (BlOb) (with the natural choice /x 2 = to 2 ) the estimate for r 3 



r 3 (TO 2 ) = |/ l f)(^ 2 ) + ^iV m (2 / 3 ) fe ^H„ 



n=0 



T(p + k + l-n) 
T(u+1- n) 



(B12) 



This gives us the values for r 3 = t 3 (to 2 ) the following estimates: 

r 3 = 1339.4(JV> = 3) , 987.3(A^/ = 4) , 690.1(JV> = 5) 



(B13) 



The principal origin of the uncertainties in these expressions is the uncertainty in the residue parameter N m (roughly 
±0.020, i.e., less than 4 %), implying an uncertainty in r 3 of a few percent (below 4 %). 

An analysis similar to this one has been performed in Ref. [72] . There, however, the term c 3 and the coefficients 

feu were not included in the analysis. 



2 the pole position in [1/1] varies 
1.00 in Nf = 5 case. 



Nonetheless, there is a possibility that at two-loop order the kinetic term contributes to an IR renormalon at b = +1 in Bg(b), 
ci. Ref. [91]. 

Though, the /^-dependence of this position is rather strong. E.g., when /t 2 varies by 10 % around m q 
between —1.6 and —0.7 in Nf = 3 case, between —1.2 and —0.7 in Nf = 4 case, and between —1.26 and 
We used for MS scheme coefficient C4 the estimated values ( B6 I, with C4 = 123.7, 97.2, 86.2, for Nf = 3, 4, 5, respectively, from Ref. 1861 . 
Simpler Pade-based estimates of C4 were obtained in Ref. 89 : C4 = 40, 70 for Nf = 4, 5, respectively (and a large negative and uncertain 
value C4 = —850 for Nf = 3). The C4 = 40 value (for Nf = 4) in this case differs substantially from the value C4 = 97.2. If we repeat 
for the C4 = 40 value [Nf = 4) the same procedure described above, we obtain C4 0.0055 (for C4 = 97.2 we got: C4 = —0.1610), hence 

5.9 (before: b a -0.96), 



the expressions of h^ 71 



(m 2 ) of Eq. (BlOb I change, and the pole of [1/1] ( an ) (6) becomes b ; 



theoretical pole b = — 1. Furthermore, from the requirement that [2/1] ( an ) (f>) have the pole at b ■■ 



-1 we now get /i.g m ^(m 2 ) 



not close to the 
4.10 



(before: —28.28), and using this value of /ig in the Pade [1/2] (an.) (&) we obtain the pole nearest to the origin b = 1.94 (before: 
b = —1.001). This indicates that the estimate C4 = 40 (for Nf = 4) is not giving results consistent with the theoretical expectations of 



the renormalon structure of B 



(an.) 
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Appendix C: Variation of pQCD coupling with scales and schemes 

In this Appendix we give the relation between eto = a p t(Qo; c 2 i c 3°^ ■ ■ ■) ana - a = a pt(Q 2 ; c 2, c 3 , . . .), where the 
latter is expressed as power expansion of the former (cf. Appendix A of Ref. [90j for details) 



a + a (—x) + a Q (x — cix + 6c 2 ) 

5 1 

+&g(— x 3 + 2°i x2 — c 2 x — Sx8c 2 + 2^ C3 ) 

4 j^ + i^i + ^ + ^x 2 

- ! ■<■■:•' - 3 Cl Sc 2 - 2Sc 3 )x + { l -c^5c 2 + ^(Sc 2 ) 2 - ^ Cl 6c 3 + ^c 4 )] + O(og) 



x 

(0) 



where we denote 



a p t(Q 2 - 1 c 2 ,c 3 ,...) , a = a pt {Ql; c 2 0) , cf ] , . . .) , 

/02 



X = Po In —2 , <5c fc = c fc - c fc 



(0) 



QV 



(Cl) 

(C2) 
(C3) 



For the purposes of our paper, it is sufficient to consider in the above relation (Cl) terms up to (including) terms 



The 3-loop threshold connection of a p t in MS scheme at the threshold scale fi 2 hv = (JCm c ) 2 (where K, ~ 1; usually 
K. = 2) can be written as the following relation between a pt (n 2 hr + 0; Nj = 4) = a + and a pt (n 2 hr — 0; Nf = 3) = a_: 



a + 



a_ [l + xia_ + x 2 a 2 _ + x 3 a 3 _ + 0(at)] , 



where 



x\ = -k\ , x 2 = -k 2 + 2k\ 



x 3 = —k 3 + hk\k 2 — bk\ 



(C4) 
(C5) 



and the coefficients fc, were calculated in Ref. 166 



ki 



-4, k 2 = —£ 2 h -—£ h + —, 
6 ' 36 71 24 72 ' 

_ 131 , 2 , (-6793 + 281 A,) 

216 h 576 '' 1728 h + 



82043 A/ „. 564731 2633 Ar 

-C(3) + - — % 



27648 



124416 31104 



(C6) 



where — \n(fi 2 hl . / 'm 2 ) = ln/C 2 , and N e in Eq. (C6| is the number of light quark flavors, i.e., N e = 3 in the considered 
case of transition from Nf = 4 to Nf = 3. 
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